home *** CD-ROM | disk | FTP | other *** search
/ Gold Fish 1 / GoldFishApril1994_CD1.img / d1xx / d199 / asimplex / matrix.c < prev    next >
Encoding:
Text File  |  1989-03-31  |  7.0 KB  |  192 lines
  1. /*****************************************************************************
  2.  * Modul                : matrix.c                                           *
  3.  * Zweck                : Funktionen für Matrixinversion, Matrizenmultipli-  *
  4.  *                        kation und Matrizenzugriff                         *
  5.  * Autor                : Stefan Förster                                     *
  6.  *                                                                           *
  7.  * Datum      | Version | Bemerkung                                          *
  8.  * -----------|---------|--------------------------------------------------- *
  9.  * 21.12.1988 |         | M()                                                *
  10.  *            |         | SetM()                                             *
  11.  *            |         | Mult()                                             *
  12.  * 05.01.1989 |         | Invers()                                           *
  13.  * 18.01.1989 |         | Mult() ohne M() und SetM()                         *
  14.  * 30.01.1989 |         | Invers() ohne Verwendung von LRSubst(),LRPart())   *
  15.  *            |         | Bug in M(): Typumwandlung (LONG) wegen Overflow    *
  16.  *            |         | Bug in SetM(): Typumwandlung (LONG) wegen Overflow *
  17.  * 06.02.1989 | 0.0     | Invers() ohne M() u. SetM(); Anpassung an ASimplex *
  18.  *            |         | Mult() angepaßt an ASimplex                        *
  19.  *            |         | M() angepaßt an ASimplex                           *
  20.  *            |         | SetM() angepaßt an ASimplex                        *
  21.  * 02.03.1989 | 0.1     | Bug in M(): if(flag==PHASE_I) statt if(PHASE_I)    *
  22.  * 07.03.1989 | 1.0     | Bug in Invers(): fabs() statt ABS() (Typ SHORT !!) *
  23.  *****************************************************************************/
  24.  
  25.  
  26.  
  27.  
  28. /*****************************************************************************
  29.  * USHORT Invers(inv,n,p)                                                    *
  30.  * Zweck     : Berechnung der inversen Matrix mit Hilfe des Gauß-Algorithmus *
  31.  * Input     : inv    - Matrix, von der die inverse Matrix berechnet werden  *
  32.  *                      soll                                                 *
  33.  *             n      - Größe der (n,n)-Matrix                               *
  34.  * Output    : inv    - inverse Matrix                                       *
  35.  *             p      - n-Permutationsvektor                                 *
  36.  * Ergebnis  : INVERTABLE/NOT_INVERTABLE                                     *
  37.  *****************************************************************************/
  38.  
  39. USHORT Invers(inv,n,p)
  40. DOUBLE *inv;
  41. SHORT  n,*p;
  42. {
  43.   REGISTER DOUBLE *ptr1, *ptr2, *ptr3, s;
  44.   REGISTER SHORT  i,j,k;
  45.   REGISTER DOUBLE max;
  46.  
  47.   for(k=1; k<=n; ++k) {
  48.     max = 0.0;
  49.     ptr2 = inv + ((LONG)(k-1)*(LONG)(n+1));
  50.     for(i=k; i<=n; ++i) {
  51.       s = 0;
  52.       ptr1 = inv + ((LONG)(i-1)*(LONG)n+(LONG)(k-1));
  53.       for(j=k; j<=n; ++j) {
  54.         s += fabs(*ptr1);
  55.         ++ptr1;
  56.       }
  57.       if(s==0.0) return(NOT_INVERTABLE);
  58.       s = fabs(*ptr2)/s;
  59.       ptr2 += n;
  60.       if(s>max) {
  61.         max = s;
  62.         p[k-1] = i;
  63.       }
  64.     }
  65.     if(max<EPS_INV) return(NOT_INVERTABLE);
  66.     if((i=p[k-1]) != k) {
  67.       ptr1 = inv+(i-1)*n;
  68.       ptr2 = inv+(k-1)*n;
  69.       for(j=0; j<n; ++j) {
  70.         s = *ptr1;
  71.         *ptr1++ = *ptr2;
  72.         *ptr2++ = s;
  73.       }
  74.     }
  75.     s = *(inv + ((LONG)(k-1)*(LONG)(n+1)));
  76.     ptr1 = inv + ((LONG)(k-1)*(LONG)n);
  77.     for(j=1; j<=n; ++j) {
  78.       if(j!=k) {
  79.         *ptr1 = -(*ptr1)/s;
  80.         ptr2 = inv + (j-1);
  81.         ptr3 = inv + (k-1);
  82.         for(i=1; i<=n; ++i) {
  83.           if(i!=k) *ptr2 += (*ptr3)*(*ptr1);
  84.           ptr2 += n;
  85.           ptr3 += n;
  86.         }
  87.       }
  88.       ++ptr1;
  89.     }
  90.     ptr1 = inv + (k-1);
  91.     for(i=1; i<=n; ++i) {
  92.       if(i!=k) *ptr1 /= s;
  93.       else     *ptr1 = 1.0/s;
  94.       ptr1 += n;
  95.     }
  96.   }
  97.  
  98.  
  99.   for(k=n-1; k>0; --k) {
  100.     if((j=p[k-1]) != k) {
  101.       ptr1 = inv + (k-1);
  102.       ptr2 = inv + (j-1);
  103.       for(i=1; i<=n; ++i) {
  104.         s = *ptr1;
  105.         *ptr1 = *ptr2;
  106.         *ptr2 = s;
  107.         ptr1 += n;
  108.         ptr2 += n;
  109.       }
  110.     }
  111.   }
  112.  
  113. return(INVERTABLE);
  114. }
  115.  
  116.  
  117.  
  118. /*****************************************************************************
  119.  * VOID Mult(A,B,erg,n,m,k)                                                  *
  120.  * Zweck     : Matrizenmultiplikation erg = A*B                              *
  121.  * Input     : A      - (n,m)-Matrix                                         *
  122.  *             B      - (m,k)-Matrix                                         *
  123.  *             n,m,k  - Zeilen- und Spaltenanzahlen                          *
  124.  * Output    : erg    - (n,k)-Matrix A*B                                     *
  125.  *****************************************************************************/
  126.  
  127. VOID    Mult(A,B,erg,n,m,k)
  128. DOUBLE  *A,*B,*erg;
  129. SHORT   n,m,k;
  130. {
  131.   REGISTER DOUBLE *Aptr,*Bptr,s;
  132.   REGISTER SHORT  h,j,i;
  133.  
  134.   for(i=0; i<n; ++i) {          /* Zeilen von erg  */
  135.     for(j=0; j<k; ++j) {        /* Spalten von erg */
  136.       Aptr = A;
  137.       Bptr = B;
  138.       for(s=0.0,h=0; h<m; ++h) {
  139.         s += (*Aptr++)*(*Bptr);
  140.         Bptr += k;
  141.       }
  142.       *erg++ = s;
  143.       ++B;
  144.     }
  145.     A += m;
  146.     B -= k;
  147.   }
  148. }
  149.  
  150.  
  151.  
  152. /*****************************************************************************
  153.  * DOUBLE M(matrix,n,i,j,flag)                                               *
  154.  * Zweck     : M() liest die Komponente a[i,j] einer Matrix                  *
  155.  * Input     : matrix - Zeiger auf die Matrix                                *
  156.  *             n      - Anzahl Spalten                                       *
  157.  *             i      - Zeilenindex                                          *
  158.  *             j      - Spaltenindex                                         *
  159.  *             flag   - PHASE_I/PHASE_II                                     *
  160.  * Bemerkung : Falls j>n und flag==PHASE_I, gibt M() Werte der angehängten   *
  161.  *             Einheitsmatrix aus (0 oder 1)                                 *
  162.  *****************************************************************************/
  163.  
  164. DOUBLE M(matrix,n,i,j,flag)
  165. DOUBLE *matrix;
  166. SHORT  n,i,j;
  167. USHORT flag;
  168. {
  169.   if(flag==PHASE_I && j>n) return(i==j-n ? 1.0 : 0.0);
  170.   else                     return(*(matrix+((LONG)(i-1)*(LONG)n+(LONG)(j-1))));
  171. }
  172.  
  173.  
  174.  
  175. /*****************************************************************************
  176.  * VOID SetM(matrix,n,i,j,value)                                             *
  177.  * Zweck     : SetM() setzt die Komponente a[i,j] einer Matrix auf 'value'   *
  178.  * Input     : matrix - Zeiger auf die Matrix                                *
  179.  *             n      - Anzahl Spalten                                       *
  180.  *             i      - Zeilenindex                                          *
  181.  *             j      - Spaltenindex                                         *
  182.  *             value  - a[i,j] = value                                       *
  183.  *****************************************************************************/
  184.  
  185. VOID    SetM(matrix,n,i,j,value)
  186. DOUBLE  *matrix;
  187. SHORT   n,i,j;
  188. DOUBLE  value;
  189. {
  190.   *(matrix+((LONG)(i-1)*(LONG)n+(LONG)(j-1))) = value;
  191. }
  192.