tools

#include <stdio.h>     //printf
#include <malloc.h>    //alloc+free
#include <stdlib.h>    //system
#include <math.h>      //fabs
#define CLR "cls"
//********************* globale Daten ***************************************
int ZA,                //****** Zeilen-Anfang
    ZE,                //****** Zeilen-Ende
    SA,                //****** Spalten-Anfang
    SE,                //****** Spalten-Ende
    *I;                //****** Index-mitLaeufer fuer X-Werte (SpaltenTausch)
double *A,             //****** Koeffizienten-Matrix
      *X,              //****** LoesungsVektor
      *Y;              //****** LoesungenVektor
//********************* Proceduren Anmeldung ********************************
void max(void);        //Suchen von Max / Tauschen der Zeilen+Spalten
void piv(void);        //rechne Fuehre PivotZeile aus
int nul(void);        //rechne Fuehre PivotZeile aus
void fhl(char);        //FehlerAusgabe
//************************ M M A A I I N N *******************************
void main(void){
int z,s,fehler,rg=0,ZM;                    //laufVar
//***************************************************************************
//********************** E I N G A B E T E I L ******************************
//***************************************************************************
system(CLR);   fhl(0);                   //****** EINGABE DER NxM
do{fhl(2);scanf("%d",&SE);if(SE<1) printf("\a");}while(SE<1);
do{fhl(1);scanf("%d",&ZE);if(ZE<1) printf("\a");}while(ZE<1);
//******************** FEHLER 0 ********************************************
if(ZE<SE){fhl(6);exit(0);}                //UnterBestimmung Fehler[1]
if(ZE*SE>6400) {fhl(3);exit(0);}          //Dimension wird zu gross
//**************************************************************************
A=(double*)malloc(ZE*SE*sizeof(double)); //Hole Matrix
X=(double*)malloc(SE* sizeof(double));   //Hole X-Vektor
I=(int *)malloc(SE* sizeof(int ));     //Hole IndexVektor
Y=(double*)malloc(ZE* sizeof(double));   //Hole LoesungenVektor
if((A==NULL)||(X==NULL)||(I==NULL)||(Y==NULL)) //beim Holen Fehler aufgetreten
    {fhl(3);exit(0);}
for(s=0;s<SE;s++) *(I+s)=s;        //Index von X-Vektor festlegen
ZA=SA=0;                           //StartWerte des Algos
ZM=ZE;                             //merke die ZeilenAnzahl
system(CLR);   fhl(4);             //***********************************
for(z=0;z<ZE;z++){                 //****** EINGABE DER MATRIX *********
for(s=0;s<SE;s++){                //***********************************
    printf("a(%d,%d)= ",z+1,s+1);   //***********************************
    scanf("%lf",(A+z*SE+s));}}      //***********************************
system(CLR);   fhl(5);             //***********************************
for(z=0;z<ZE;z++){                 //*** EINGABE DES LÖSUNGENVEKTORS ***
printf("y(%d)= ",z+1);            //***********************************
scanf("%lf",(Y+z));}              //***********************************
//********************** Rechen-Schleife ************************************
do{
   fehler=nul();                    //Pruefe auf Null-Zeilen
   if(fehler!=0) continue;          //...
   max();                           //suche Max-PivoElement
   //printf("nach Max"); for(z=0;z<ZE;z++){printf("\n");for(s=0;s<SE;s++){
   //printf("%le ",*(A+(z*SE+s)));} printf("%le ",*(Y+z));} getchar();
   piv();                           //eliminiere Spalte
   rg++;                            //Rang der Matrix jetzt
//   printf("nach elim");for(z=0;z<ZE;z++){printf("\n");for(s=0;s<SE;s++){
//   printf("%le ",*(A+(z*SE+s)));} printf("%le ",*(Y+z));} getchar();
   ZA++;SA++;                       //Matrix wird kleiner
}while((ZA<ZE)&&(SA<SE));          //Matrix Pruefen ob Ende erreicht
//*************************** Fehler 1 *************************************
if(fehler==2) {fhl(7);exit(0);}
if(fehler==0) {fehler=nul(); }//erkenne die letzte Zeile
if(fehler==2) {fhl(7);exit(0);}
if(rg<SE)     {fhl(8);exit(0);}
//******************************* Errechnen des X-Vektors *******************
for(s=SE-1;s>=0;s--){
*(X+s)=*(Y+s)/(*(A+SE*s+s));       //bestimme das X-Element [X = Y div Koeff]
for(z=SE-2;z>=0;z--){              //von allen restlichen Y's wird Koeff*X
    *(Y+z)-=(*(A+SE*z+s))*(*(X+s));}}//gleich abgezogen
//*************************** Ausgabe ***************************************
system(CLR);

if(ZE<ZM) printf("\n\t\tGleichungsSystem war ueberbestimmt.\n");
if(ZE==ZM) printf("\n\t\tGleichungsSystem war optimal bestimmt.\n");
for(s=0;s<SE;s++){                //ueber alle X's
    for(z=0;*(I+z)!=s;z++);        //suche Index <S> == alte Ordnung
    if(*(X+z)<0.0) printf("x(%d) = %.5lf \t== %.15le\n",s+1,*(X+z),*(X+z)); //Ausgabe
    if(*(X+z)>=0.0) printf("x(%d) = %.5lf \t== %.15le\n",s+1,*(X+z),*(X+z));} //Ausgabe
free(A); free(X); free(I); free(Y);} //Freigabe der 'Mallox'
//********************** E N D E   M A I N **********************************

//***************************************************************************
//** UP1 ********* Suchen MAX: und Tauschen der Spalten und Zeilen *********
//***************************************************************************
void max(void){          //Suchen von Max / Tauschen der Zeilen+Spalten
int s,z,sm=SA,zm=ZA;     //default
double tmp=*(A+ZA*SE+SA);//default
for(z=ZA;z<ZE;z++){
   for(s=SA;s<SE;s++){
     if(fabs(*(A+z*SE+s))>fabs(tmp)) {sm=s;zm=z;tmp=*(A+z*SE+s);}}}
for(s=SA;s<SE;s++){     //**** Zeilen-Tausch Matrix A
    tmp=*(A+ZA*SE+s);    //merke das ZeilenElement aus 1.Zeile
    *(A+ZA*SE+s)=*(A+zm*SE+s); //Überschreibe das ZeilenElement 1.Zl.mit MAX
    *(A+zm*SE+s)=tmp;}   //MAX-Zeile-Elem.= mit 1.Zeile
tmp=*(Y+ZA);          //**** Zeilen-Tausch Vektor Y
*(Y+ZA)=*(Y+zm);      //Überschreibe das ZeilenElement 1.Zl.mit MAX
*(Y+zm)=tmp;          //MAX-Zeile-Elem.= mit 1.Zeile
for(z=0;z<ZE;z++){    //**** Spalten-Tausch alle SpaltenElemente
    tmp=*(A+z*SE+SA); //merke das SpaltenElement aus 1.Spalte
    *(A+z*SE+SA)=*(A+z*SE+sm); //Überschreibe das SpaltenElement 1.Sp.mit MAX
    *(A+z*SE+sm)=tmp;} //MAX-Spalte-Elem.= mit 1.Spalte
zm=*(I+SA);           //**** Spalten-Tausch Vektor- I (X-Index)
*(I+SA)=*(I+sm);      //nehme zm als tmp, da *I == int; brauch zm nicht mehr
*(I+sm)=zm;}           //**********************

//***************************************************************************
//** UP2 ********* Fuehre die PivotZeile aus; mit Probe auf 0.0 **************
//***************************************************************************
void piv(void){
int s,z,l1,l2,l3;                   //laufVar
double fak;                         //Pivot-Faktor
for(z=ZA+1;z<ZE;z++){              //ueber alle Zeilen
   fak=*(A+z*SE+SA)/(*(A+ZA*SE+SA));//1.Sp-Elem div Pivot
   l1=l2=l3=0;                      //Exponent setzen; auf NUL Merken
   if(*(Y+z)!=0.0) {                //*** Exponent Minuend
     l1=(int)fabs(log10(fabs((*(Y+z)))));}
   if(fak*(*(Y+ZA))!=0.0) {         //*** Exponent Subtrahend
     l2=(int)fabs(log10(fabs(fak*(*(Y+ZA)))));}
   *(Y+z) -=fak*(*(Y+ZA));          //Y-Vektor Anpassen
   if(*(Y+z)!=0.0) {                //*** Exponent Ergebnis
     l3=(int)fabs(log10(fabs((*(Y+z)))));}
   if( (((l1-l2)*(l1-l2))<2) &&     //*** bei starker Exponent-
         (((l1-l3)*(l1-l3))>25)) *(Y+z)=0.0;    //*** Abweichung NULL setzen
   for(s=SA;s<SE;s++){               //ueber alle Spalten
     l1=l2=l3=0;                     //Exponent setzen auf NUL Merken
     if(*(A+z*SE+s)!=0.0) {          //*** Exponent Minuend
       l1=(int)fabs(log10(fabs((*(A+z*SE+s)))));}
     if(fak*(*(A+(ZA*SE+s)))!=0.0) { //*** Exponent Subtrahend
       l2=(int)fabs(log10(fabs(fak*(*(A+(ZA*SE+s))))));}
     //****** printf("%le %le %d %d",*(A+z*SE+s),fak*(*(A+(ZA*SE+s))),l1,l2);
     *(A+z*SE+s)-=fak*(*(A+(ZA*SE+s))); //Ausfuehren des Pivotzeile
     if(*(A+z*SE+s)!=0.0) {             //*** Exponent Ergebnis
       l3=(int)fabs(log10(fabs((*(A+z*SE+s)))));}
     //****** printf(" %le %d\n",*(A+z*SE+s),l3);
     if( (((l1-l2)*(l1-l2))<2) &&                  //*** bei starker Exponent
         (((l1-l3)*(l1-l3))>25)) *(A+z*SE+s)=0.0;}}//*** Abweichung NULL setzen
} //Spalten abgleichen
//***************************************************************************
//** UP3 ********* Loeschen der Nul-Zeilen ***********************************
//***************************************************************************
int nul(void){
int s,z,ma,my,ret=0;                 //LaufVar
double tmp;                          //AblegeVar
for(z=ZA;z<ZE;z++){                 //Über alle Zeilen
   ma=0;my=0;                        //Zaehl-Var loeschen
   for(s=SA;s<SE;s++){               //zaehle die Nullen in der Zeile
     if(*(A+z*SE+s)==0.0) ma++;}     //wenn gefunden dann Zaehle hoch
   if(*(Y+z)==0.0) my++;             //ist Y==0 oder nich?
   if((ma==(SE-SA))&&(my==0)){       //Unterbestimmung; grober Fehler
     SA=SE;return(2);}               //bereite das Ende vor und Tschuess
   if((ma==(SE-SA))&&(my==1)){       //sind alle Elemente der Zeile == 0
     ZE--;                           //ZeilenEnde neu anpassen
     for(s=SA;s<SE;s++){        //**** Zeilen-Tausch der Matrix A
      tmp=*(A+z*SE+s);          //merke das ZeilenElement aus aktuelle Zeile
      *(A+z*SE+s)=*(A+ZE*SE+s); //Überschreibe das ZeilenElement
      *(A+ZE*SE+s)=tmp;}        //MAX-Zeile-Elem.= mit 1.Zeile
     tmp=*(Y+z);                //**** Zeilen-Tausch Vektor Y
     *(Y+z)=*(Y+ZE);            //Überschreibe das ZeilenElement
     *(Y+ZE)=tmp;               //....
     z--;ret=1;}}               //aktuelle Zeile ist getauscht nochMal Zeile
return(ret);}                   //alles ok return(0);or(1);
//***************************************************************************
//** UP4 ********* Ausgabe von Saetzen und Fehlern ***************************
//***************************************************************************
void fhl(char nummer){
char *S[9]={           //Text-Ausgaben
"                            Gauss-Algo\n"                            //0
"                           ===========\n"                           //0
"                            _ _\n"                                 //0
"                        A * X= Y\n"                                 //0
"geg. sind Koeffizienten Matrix A (mit z-Zeilen & s-Spalten) und der\n"
"Loesungs-Vektor Y (mit z-Zeilen). Gesucht ist der Vektor X (mit s-Zeilen),\n"
"der dieses Gleichungssystem loest.\n\n"                              //0
"    Ú x[1]*a[1,1]   x[2]*a[1,2] ... x[s]*a[1,s] ¿ = Y[1]\n"       //0
"    ³ x[1]*a[2,1]   x[2]*a[2,2] ... x[s]*a[2,s] ³ = Y[2]\n"       //0
"    : :      :      :     :     ...   :    :     : = : \n"       //0
"    À x[1]*a[z,1]   x[2]*a[z,2] ... x[s]*a[z,s] Ù = Y[z]\n\n",    //0
"\nGib die ZeilenZahl von A z:= ",                                 //1
"\nGib die SpaltenZahl von A s:= ",                                 //2
"\a Matrix ist zu gross, zu wenig Speicher <EXIT>\n",                 //3
"\n gib die Werte der MATRIX A ein\n"                                //3
" Elemente-Format a(z,s), d.h. a(zeile,spalte)\n\n",                 //4
"\n gib die Werte des LoesungsVektors Y ein\n"                        //5
" Elemente-Format y(z), d.h. y(zeile)\n\n",                          //5
"\n\n\a\tDas GleichungsSystem ist unterbestimmt, zur korrekten Aussage\n"
      "\tueber den Vektor- X muessen mindestens genausoviele Zeilen wie\n"
      "\tSpalten vorhanden sein. Das GleichungsSystem besitzt unendlich\n"
      "\tviele Loesungen. <zu wenig Zeilen angegeben>\n",            //6
"\n\n\a\tDas Gleichungs- System kann nicht geloest werden. Die um\n"//7
      "\tY- Vektor erweiterte Koeffizienten- Matrix A hat einen\n"//7
      "\thoeheren Rang, als nur die Matrix A, dh. rg(A,Y) > rg(A),\n"//7
      "\tueber ein und den selben Vektor werden zwei oder mehr ver-\n"//7
      "\tschiedene Aussagen gemacht. Fuer z=2 und s=1 ein Bsp.\n"//7
      "\t2*x=2 und 3*x=1\n",                                         //7
"\n\n\a\tDas GleichungsSystem kann nicht eindeutig geloest werden.\n" //8
      "\tEs gibt unendlich viele Loesungen. Die Gleichungen sind\n"   //8
      "\tlinear abhaengig.\n"};                                       //8
printf("%s",S[nummer]);}
//***************************************************************************
//********************** ENDE PROGRAM GAUSS - ALGO **************************
//***************************************************************************
//***************************************************************************
//****** p r o g r a m m i e r t   v o n   d a v e   s u n   i n   14h ******
//***************************************************************************
//***************************************************************************
//************************** 30.11.1994 by Dave Sun *************************
//***************************************************************************