/********************************************************************/
 /****                                                            ****/
 /****                                                            ****/
 /****    Program          : MatrixI.c                            ****/
 /****                                                            ****/
 /****    Version          : 01.00                                ****/
 /****                                                            ****/
 /****    Erstversion      : 01.07.1990                           ****/
 /****                                                            ****/
 /****    Letzte Änderung  : 08.08.1990                           ****/
 /****                                                            ****/
 /****    Compiliert mit   : siehe MAKEFILE                       ****/
 /****                                                            ****/
 /****    Gelinkt mit      : siehe MAKEFILE                       ****/
 /****                                                            ****/
 /********************************************************************/
 /****                                                            ****/
 /****                                                            ****/
 /****               Copyright by Rüdiger Dreier                  ****/
 /****                                                            ****/
 /****                                                            ****/
 /********************************************************************/
 
 /* Einbinden der Include-Files */
 #include <exec/types.h>
 #include <string.h>
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
 
 LONG Zeilen,MZeilen;
 LONG Spalten,MSpalten;
 
 DOUBLE Zeile[20];
 DOUBLE MatrixE[20][20];
 DOUBLE MatrixL[20][20];
 DOUBLE MatrixA[20][20];
 char NULLSpalte[20];
 
 SHORT Loesen(VOID)
  {
   SHORT i,l,Loesbarkeit,Ziel=0;
   
   MZeilen=Zeilen;
   MSpalten=0;
   
   /* Erster Schritt: Alle NullSpalten eliminieren und merken */
   for(i=0;i<Spalten;i++)
    {
     if(!NullSpalte(i))
      {
       NULLSpalte[i]=0;
       for(l=0;l<Zeilen;l++)
        {
         MatrixL[l][MSpalten]=MatrixE[l][i];
        }
       MSpalten++;
      }
     else
      {
       NULLSpalte[i]=1;
      }
    }
   
   /* Zweiter Schritt: Lösen */
   Loesbarkeit=MLoesen();
   
   for(i=0;i<Spalten;i++)
    {
     for(l=0;l<Zeilen;l++)
      {
       MatrixA[l][i]=0.0; /* Löschen */
      }
    }
   /* Dritter Schritt: NullSpalten wieder herstellen */
   for(i=0;i<MSpalten;i++)
    {
     if(NULLSpalte[i])
      {
       Ziel++;
      }
     for(l=0;l<Zeilen;l++)
      {
       MatrixA[l][Ziel]=MatrixL[l][i];
      }
     Ziel++;
    }
   return(Loesbarkeit);
  }
 
 SHORT MLoesen(VOID)
  {
   SHORT k,flag=0,i,akt,mn;
   DOUBLE HILFE,HILFE2;
   
   mn=min(MZeilen,MSpalten-1)-1;
   for(k=0;k<mn;k++)
    {
     if(MatrixL[k][k]==NULL)flag=Tauschen(k);
     
     if(MatrixL[k][k]!=NULL)
      {
       /* Alle weiteren MZeilen in der aktuellen Spalte auf NULL bringen */
       HILFE =MatrixL[k][k];
       for(akt=0;akt<MZeilen;akt++)
        {
         if(akt!=k)
          {
           HILFE2=MatrixL[akt][k];
           for(i=0;i<MSpalten;i++)
            {
             MatrixL[akt][i]=MatrixL[akt][i]*HILFE-MatrixL[k][i]*HILFE2;
            }
          }
        }
      }
    }
   if(MatrixL[k][k]==NULL)Tauschen(k);
   
   if(MatrixL[k][k]!=NULL)
    {
     HILFE =MatrixL[k][k];
     for(akt=0;akt<MZeilen;akt++)
      {
       if(akt!=k)
        {
         HILFE2=MatrixL[akt][k];
         for(i=0;i<MSpalten;i++)
          {
           MatrixL[akt][i]=MatrixL[akt][i]*HILFE-MatrixL[k][i]*HILFE2;
          }
        }
      }
    }
   /* Diagonale auf 1 bringen */
   for(k=0;k<MZeilen;k++)
    {
     if(MatrixL[k][k]!=NULL)
      {
       for(i=k+1;i<MSpalten;i++)
        {
         MatrixL[k][i]/=MatrixL[k][k];
        }
       MatrixL[k][k]=1.00;
      }
    }
   /* Die Gleichung ohne NullSpalten ist evt. eindeutig lösbar, */
   /* die mit Nullspalten aber auf gar keinen Fall              */
   if(MZeilen>=MSpalten-1)
    {
     flag=NULL;
    }
   else 
    {
     flag=2; /* Nicht eindeutig lösbar */
    }
   if(MSpalten<Spalten)flag=2; /* Es gab NullSpalten */
   
   mn=min(MZeilen+1,MSpalten)-1;
   for(k=0;k<mn;k++)
    {
     SHORT flag2;
     flag2=1;
     if(MatrixL[k][k]==NULL && MatrixL[k][MSpalten-1]!=NULL)
      {
       for(i=k+1;i<MSpalten-1;i++)
        {
         /* Das Diagonalelement kann zwar NULL sein, es kann aber sein, daß */
         /* noch ein weiteres Element != NULL vorhanden ist.                */
         if(MatrixL[k][i]!=NULL)
          {
           flag2=2;
          }
        }
       flag=flag2; /* DiagonalElement der Matrix 0, Element des LösungsVektor !=0 */
       if(flag==1)break; /* Nicht mehr lösbar */
      }
     if(MatrixL[k][k]==NULL && MatrixL[k][MSpalten-1]==NULL)
      {
       flag= 2; /* Beide Elemente = 0 */
      }
    }
   for(k=mn;k<MZeilen;k++)
    {
     /* Alle Zeilen, die nicht mehr in der Diagonalen enthalten sind, */
     /* müssen ganz NULL sein, sonst nicht lösbar.                    */
     if(MatrixL[k][MSpalten-2]==NULL && MatrixL[k][MSpalten-1]!=NULL)
      {
       flag= 1;
      }
    }
   return(flag);
  }
 
 
 SHORT Tauschen(SHORT k)
  {
   SHORT flag,l,i;
   DOUBLE HILFE;
   /* Tauschen */
   flag=1;
   
   for(l=k+1;l<MZeilen;l++)
    {
     if(MatrixL[l][k]!=NULL)
      {
       for(i=0;i<MSpalten;i++)
        {
         HILFE=MatrixL[l][i];
         MatrixL[l][i]=MatrixL[k][i];
         MatrixL[k][i]=HILFE;
         flag=0;
        }
       l=MZeilen;
      }
    }
   return(flag);
  }
 
 SHORT NullSpalte(SHORT k)
  {
   SHORT i,fl=1;
   for(i=0;i<Zeilen;i++)
    {
     if(MatrixE[i][k]!=0.0)fl=0;
    }
   return(fl);
  }