Lezione del 10/7/2006

  • Esempio di un programma complesso, soluzione di un sistema lineare tramite la fattorizzazione LU di una matrice quadrata

    File ex025.c

      1/*
  2 * Esempio "complesso", fattorizzazione LU di Gauss di una Matrice
  3 */
  4
  5#include <stdio.h>
  6#include <string.h>
  7#include <stdlib.h>
  8#include <math.h>
  9
 10/* prototipo della funzione */
 11int fattorizzazione_LU( double [], int, int, int [] ) ;
 12
 13
 14#define A(I,J) A[(I)+(J)*lda]
 15
 16int /* valore di ritorno se 0 = tutto ok */
 17    /* se 1 = matrice singolare */
 18fattorizzazione_LU( double A[],       /* matrice da fattorizzare memorizzata per colonne */
 19                    int    N,         /* dimensione della matrice da fattorizzare */
 20                    int    lda,       /* numero di righe usate in allocazione     */
 21                    int    perm[] ) { /* vettori di interi di domensione >= N che
 22                                         conterra' la permutazione */
 23  int    i, j, k, imax ;
 24  double amax, acheck ;
 25
 26  /* Inizializzo il vettore della permutazione */
 27  for ( i = 0 ; i < N ; ++i ) perm[i] = i ;
 28
 29  for ( i = 0 ; i < N-1 ; ++i ) {
 30    /* ricerca del pivot, cerco elemento di modulo massimo 
 31       nella colonna i-esima, cioe' A(i,i), A(i+1,i), ..., A(N-1,i) */
 32    imax = i ;
 33    amax = fabs( A(i,i) ) ;
 34    for ( j=i+1 ; j < N ; j++ ) {
 35      acheck = fabs(A(j,i));
 36      if ( acheck > amax ) {
 37        acheck = amax ;
 38        imax   = j ;
 39      }
 40    }
 41
 42    /* controllo che acheck > 0, se 0 la matrice e' singolare */
 43    if ( acheck == 0 ) return 1 ;
 44    if ( imax != i ) { /* scambio la riga i con la riga imax */
 45      /* scambio indici sul vettore della permutazione */
 46      k          = perm[i] ;
 47      perm[i]    = perm[imax] ;
 48      perm[imax] = k ;
 49      /* scambio delle righe nella matrice A */
 50      for ( j=0 ; j < N ; ++j ) {
 51        amax      = A(i,j) ;
 52        A(i,j)    = A(imax,j) ;
 53        A(imax,j) = amax ;
 54      }
 55    }
 56
 57    /* costruisco la parte della matrice Li e la applico alla matrice
 58       LU in costruzione */
 59    for ( k = i+1 ; k < N ; ++k ) {
 60      amax   = A(k,i) / A(i,i) ;
 61      A(k,i) = amax ;
 62      for ( j= i+1 ; j < N ; ++j )
 63        A(k,j) -= amax * A(i,j) ;
 64    }
 65  }
 66  return 0 ; /* tutto ok */
 67}
 68
 69void
 70stampa_L( FILE * fd,
 71          double A[],    /* matrice da fattorizzare memorizzata per colonne */
 72          int    N,      /* dimensione della matrice da fattorizzare */
 73          int    lda ) { /* numero di righe usate in allocazione     */
 74  int i, j ;
 75  for ( i = 0 ; i < N ; ++i ) {	
 76    for ( j = 0 ; j < i ; ++j ) {
 77      fprintf( fd, "%-10lg ",A(i,j)) ;
 78    }
 79    fprintf( fd, "1\n") ; /* uso 10 caratteri per numero */
 80  }
 81}
 82
 83void
 84stampa_U( FILE * fd,
 85          double A[],    /* matrice da fattorizzare memorizzata per colonne */
 86          int    N,      /* dimensione della matrice da fattorizzare */
 87          int    lda ) { /* numero di righe usate in allocazione     */
 88  int i, j ;
 89  for ( i = 0 ; i < N ; ++i ) {	
 90    for ( j = 0 ; j < i ; ++j ) {
 91      fprintf( fd, "           ") ; /* 11 spazi bianchi */
 92    }
 93    for ( j = i ; j < N ; ++j ) {
 94      fprintf( fd, "%-10lg ",A(i,j)) ;
 95    }
 96    fprintf( fd, "\n") ; /* vado a capo a fine riga */
 97  }
 98}
 99
100
101#undef A /* elimino la definizione di A(I,I) */
102
103#define L(I,J) LU[(I)+(J)*lda]
104#define U(I,J) LU[(I)+(J)*lda]
105
106void
107solve_LU( double LU[],      /* matrice fattorizzata memorizzata per colonne */
108          int    N,         /* dimensione della matrice da fattorizzare */
109          int    lda,       /* numero di righe usate in allocazione     */
110          int    perm[],    /* vettori di interi di domensione >= N che
111                               conterra' la permutazione */
112          double b[],       /* temine noto */
113          double x[] ) {    /* soluzione */
114  
115  int i, j ;
116  /* applico la permutazione */
117  for ( i = 0 ; i < N ; ++i ) x[perm[i]] = b[i] ;
118
119  /* Risolvo Lz = Pb */
120  for ( i = 1 ; i < N ; ++i )
121   for ( j = 0 ; j < i ; ++j )
122     x[i] -= L(i,j)*x[j];
123
124  /* Risolvo Ux = Lz */
125  i=N;
126  do {
127   --i ;
128   for ( j = i+1 ; j < N ; ++j )
129     x[i] -= U(i,j)*x[j];
130   x[i] /= U(i,i) ; 
131  } while ( i > 0 ) ;
132}
133
134int
135main() {
136
137  double M[] = {
138    4, 3, 2, 1, /* prima colonna   */
139    1, 2, 3, 4, /* seconda colonna */
140    1, 0, 0, 1, /* terza colonna   */
141    0, 1, 2, 5  /* quarta colonna  */
142  } ;
143
144  double b[] = { 9, 11, 16, 32 } ; /* termine noto se la solzione e' 1,2,3,4 */
145  double x[4] ; /* vettore che conterra' la soluzione */
146
147  int i, perm[4] ;
148
149  int iflag, dim ;
150
151  dim = 4 ;
152
153  iflag = fattorizzazione_LU( M, dim, dim, perm ) ;
154
155  printf( "Valore del flag=%d\n",iflag);
156
157  #define A(I,J) M[(I)+(J)*4]
158  
159  /* stampa della matrice L ed U e la permutazione */
160  printf("MATRICE L\n") ;
161  stampa_L( stdout, M, dim, dim );
162
163  printf("MATRICE U\n") ;
164  stampa_U( stdout, M, dim, dim );
165
166  printf("PERMUTAZIONE\n") ;
167  for ( i = 0 ; i < dim ; ++i )
168    printf( "perm[%d]=%d\n", i, perm[i]) ;
169
170  solve_LU( M, dim, dim, perm, b, x ) ;
171
172  printf("SOLUZIONE\n") ;
173  for ( i = 0 ; i < dim ; ++i )
174    printf( "x[%d]=%lg\n", i, x[i]) ;
175
176  return 0 ;
177
178}

  • Esempio di un programma complesso, stesso problema ma soluzione cin l’uso della libreria LAPACK

    File ex026.c

      1/*
  2 * Esempio "complesso", fattorizzazione LU di Gauss di una Matrice
  3 * Ma usando le LAPACK/BLAS
  4 *
  5 */
  6
  7#include <stdio.h>
  8#include <string.h>
  9#include <stdlib.h>
 10#include <math.h>
 11
 12/*
 13   FATTORIZZAZIONE LU, PROTOTIPO FORTRAN DELLE LAPACK
 14
 15 NAME
 16      DGETRF - compute an LU factorization of a general M-by-N
 17      matrix A using partial pivoting with row interchanges
 18
 19 SYNOPSIS
 20      SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
 21
 22          INTEGER        INFO, LDA, M, N
 23
 24          INTEGER        IPIV( * )
 25
 26          DOUBLE         PRECISION A( LDA, * )
 27
 28 PURPOSE
 29      DGETRF computes an LU factorization of a general M-by-N
 30      matrix A using partial pivoting with row interchanges.
 31
 32      The factorization has the form
 33         A = P * L * U
 34      where P is a permutation matrix, L is lower triangular with
 35      unit diagonal elements (lower trapezoidal if m > n), and U
 36      is upper triangular (upper trapezoidal if m < n).
 37
 38      This is the right-looking Level 3 BLAS version of the algo-
 39      rithm.
 40
 41 ARGUMENTS
 42      M       (input) INTEGER
 43              The number of rows of the matrix A.  M >= 0.
 44
 45      N       (input) INTEGER
 46              The number of columns of the matrix A.  N >= 0.
 47
 48      A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 49              On entry, the M-by-N matrix to be factored.  On
 50              exit, the factors L and U from the factorization A =
 51              P*L*U; the unit diagonal elements of L are not
 52              stored.
 53
 54      LDA     (input) INTEGER
 55              The leading dimension of the array A.  LDA >=
 56              max(1,M).
 57
 58      IPIV    (output) INTEGER array, dimension (min(M,N))
 59              The pivot indices; for 1 <= i <= min(M,N), row i of
 60              the matrix was interchanged with row IPIV(i).
 61
 62      INFO    (output) INTEGER
 63              = 0:  successful exit
 64              < 0:  if INFO = -i, the i-th argument had an illegal
 65              value
 66
 67              > 0:  if INFO = i, U(i,i) is exactly zero. The fac-
 68              torization has been completed, but the factor U is
 69              exactly singular, and division by zero will occur if
 70              it is used to solve a system of equations.
 71
 72*/
 73
 74/* 
 75  devo definire il prototipo da chiamare in C
 76
 77
 78 */
 79
 80
 81#define F77NAME(A) A##_
 82
 83
 84extern void F77NAME(dgetrf)( int    * M, /* gli argomenti sono sempre passati per indirizzo */
 85                             int    * N,
 86                             double * A,
 87                             int    * LDA,
 88                             int    * IPIV,
 89                             int    * INFO ) ;
 90
 91
 92
 93
 94/*
 95
 96 NAME
 97      DGETRS - solve a system of linear equations  A * X = B or A'
 98      * X = B with a general N-by-N matrix A using the LU factori-
 99      zation computed by DGETRF
100
101 SYNOPSIS
102      SUBROUTINE DGETRS( TRANS, N, NRHS, A, LDA, IPIV, B, LDB,
103                         INFO )
104
105          CHARACTER      TRANS
106
107          INTEGER        INFO, LDA, LDB, N, NRHS
108
109          INTEGER        IPIV( * )
110
111          DOUBLE         PRECISION A( LDA, * ), B( LDB, * )
112
113 PURPOSE
114      DGETRS solves a system of linear equations
115         A * X = B  or  A' * X = B with a general N-by-N matrix A
116      using the LU factorization computed by DGETRF.
117
118 ARGUMENTS
119      TRANS   (input) CHARACTER*1
120              Specifies the form of the system of equations:
121              = 'N':  A * X = B  (No transpose)
122              = 'T':  A'* X = B  (Transpose)
123              = 'C':  A'* X = B  (Conjugate transpose = Transpose)
124
125      N       (input) INTEGER
126              The order of the matrix A.  N >= 0.
127
128      NRHS    (input) INTEGER
129              The number of right hand sides, i.e., the number of
130              columns of the matrix B.  NRHS >= 0.
131
132      A       (input) DOUBLE PRECISION array, dimension (LDA,N)
133              The factors L and U from the factorization A = P*L*U
134              as computed by DGETRF.
135
136      LDA     (input) INTEGER
137              The leading dimension of the array A.  LDA >=
138              max(1,N).
139
140      IPIV    (input) INTEGER array, dimension (N)
141              The pivot indices from DGETRF; for 1<=i<=N, row i of
142              the matrix was interchanged with row IPIV(i).
143
144 (LDB,NRHS)
145      B       (input/output) DOUBLE PRECISION array, dimension
146              On entry, the right hand side matrix B.  On exit,
147
148              the solution matrix X.
149
150      LDB     (input) INTEGER
151              The leading dimension of the array B.  LDB >=
152              max(1,N).
153
154      INFO    (output) INTEGER
155              = 0:  successful exit
156              < 0:  if INFO = -i, the i-th argument had an illegal
157              value
158
159*/
160
161extern void F77NAME(dgetrs)( char   * TRANS, /* gli argomenti sono sempre passati per indirizzo */
162                             int    * N,
163                             int    * NRHS,
164                             double * A,
165                             int    * LDA,
166                             int    * IPIV,
167                             double * B,
168                             int    * LDB,
169                             int    * INFO ) ;
170
171#define A(I,J) A[(I)+(J)*lda]
172
173void
174stampa_L( FILE * fd,
175          double A[],    /* matrice da fattorizzare memorizzata per colonne */
176          int    N,      /* dimensione della matrice da fattorizzare */
177          int    lda ) { /* numero di righe usate in allocazione     */
178  int i, j ;
179  for ( i = 0 ; i < N ; ++i ) {	
180    for ( j = 0 ; j < i ; ++j ) {
181      fprintf( fd, "%-10lg ",A(i,j)) ;
182    }
183    fprintf( fd, "1\n") ; /* uso 10 caratteri per numero */
184  }
185}
186
187void
188stampa_U( FILE * fd,
189          double A[],    /* matrice da fattorizzare memorizzata per colonne */
190          int    N,      /* dimensione della matrice da fattorizzare */
191          int    lda ) { /* numero di righe usate in allocazione     */
192  int i, j ;
193  for ( i = 0 ; i < N ; ++i ) {	
194    for ( j = 0 ; j < i ; ++j ) {
195      fprintf( fd, "           ") ; /* 11 spazi bianchi */
196    }
197    for ( j = i ; j < N ; ++j ) {
198      fprintf( fd, "%-10lg ",A(i,j)) ;
199    }
200    fprintf( fd, "\n") ; /* vado a capo a fine riga */
201  }
202}
203
204#undef A /* elimino la definizione di A(I,I) */
205
206int
207main() {
208
209  double M[] = {
210    4, 3, 2, 1, /* prima colonna   */
211    1, 2, 3, 4, /* seconda colonna */
212    1, 0, 0, 1, /* terza colonna   */
213    0, 1, 2, 5  /* quarta colonna  */
214  } ;
215
216  double b[] = { 9, 11, 16, 32 } ; /* termine noto se la solzione e' 1,2,3,4 */
217  double x[4] ; /* vettore che conterra' la soluzione */
218
219  int i, perm[4], one ;
220
221  int iflag, dim ;
222
223  dim = 4 ;
224  one = 1 ;
225  
226  /* richiamo il fattorizzatore delle LAPACK */
227  F77NAME(dgetrf)( &dim, &dim, M, &dim, perm, &iflag );
228
229  printf( "Valore del flag=%d\n",iflag);
230
231  #define A(I,J) M[(I)+(J)*4]
232  
233  /* stampa della matrice L ed U e la permutazione */
234  printf("MATRICE L\n") ;
235  stampa_L( stdout, M, dim, dim );
236
237  printf("MATRICE U\n") ;
238  stampa_U( stdout, M, dim, dim );
239
240  printf("PERMUTAZIONE\n") ;
241  for ( i = 0 ; i < dim ; ++i )
242    printf( "perm[%d]=%d\n", i, perm[i]) ;
243
244  for ( i = 0 ; i < dim ; ++i ) x[i] = b[i] ;
245  
246  /* richiamo il solutore delle LAPACK */
247  F77NAME(dgetrs)( "N", &dim, &one, M, &dim, perm, x, &dim, &iflag );
248
249  printf( "Valore del flag=%d\n",iflag);
250
251  printf("SOLUZIONE\n") ;
252  for ( i = 0 ; i < dim ; ++i )
253    printf( "x[%d]=%lg\n", i, x[i]) ;
254
255  return 0 ;
256
257}
258
259/*
260 * Su una macchimna UNIX si comile con
261 *
262 * gcc -Wall ex026.c -o ex026 -llapack -lblas -lm
263 *
264 * Su un mac con i Framework
265 *
266 * gcc -Wall ex026.c -o ex026 -framework Accelerate
267 *
268 */