/*     Computing formal-algebraic solution to interval linear systems
         in Kaucher arithmetic  by single-step stationary iterative     
	   method  with  modified  real  splitting  of the matrix,
		      (c) Sergey P. Shary, 1995                             */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <conio.h>
#include <string.h>
#include <dos.h>


main(int argc, char *argv[])

{   char   DataFile[64], *Z;
     int   ier, *ip, i, j, k, m, ni, Dim, DIM;
  double   g0, g1, h0, h1, p, q, s0, s1, t0, t1, u0, u1, v0, v1, Eps, Tau;
  double   *G, *d, *H, *x, *xx;  
    FILE   *fi, *fopen();
    long   IterLim;

textmode(C80); clrscr(); textcolor(RED); cprintf("\r\n\n     *** ");
textcolor(LIGHTRED);
cprintf("  Computing  formal  solutions  to  interval  linear  systems  ");
textcolor(RED); cprintf(" *** \r\n      ** "); textcolor(LIGHTRED);
cprintf("  by iterative method with real splitting, (c) S.Shary, 1995  ");
textcolor(RED); cprintf(" ** \n"); delay(1200); textcolor(LIGHTGRAY);

if(argc == 1)
{ cprintf("\n\rProgram RE_SPLIT  computes  formal  solutions  to interval");
  cprintf(" linear systems  in\n");
  cprintf("\r Kaucher arithmetic  and implements stationary one-step");
  cprintf(" iterative method  based\n");
  cprintf("\r on  modified  real splitting  of the matrix  of the system,");
  cprintf(" which  is supposed\n");
  cprintf("\r to be square. \n\n");
  cprintf("\rUsage:"); textcolor(WHITE); 
  cprintf("   RE_SPLIT  DataFile \n\n "); textcolor(LIGHTGRAY);
  cprintf("\rIn the DataFile, the following parameters  must occupy  the"); 
  cprintf(" first four positions (in order): \n");
  textcolor(WHITE); cprintf("\r %15s"," * "); textcolor(LIGHTGRAY); 
  cprintf("dimension of the interval equation (a positive integer), \n");
  textcolor(WHITE); cprintf("\r %15s"," * "); textcolor(LIGHTGRAY);
  cprintf("relative defect of the approximate solution (a positive real),\n");
  textcolor(WHITE); cprintf("\r %15s"," * "); textcolor(LIGHTGRAY);
  cprintf("an auxiliary real parameter (may be arbitrary),\n");
  textcolor(WHITE); cprintf("\r %15s"," * "); textcolor(LIGHTGRAY);
  cprintf("restriction to the number of iterations (a positive integer). \n\n");
  cprintf("\rNext the interval matrix (along the rows) and");
  cprintf(" right-hand side vector are written out in the DataFile.\n"); 
  goto final;

}
else
  strcpy(DataFile,argv[1]);

  /*  reading input data,
		   forming the "middle" matrix G for the given system      */

                          fi=fopen(DataFile,"r");
if( fi==NULL )
{  printf("\n\r    Error reading the data file %10s ! \n",DataFile);
   sound(120); delay(1000); nosound(); goto final;
}

if( fscanf(fi,"%d %lf %lf %lu",&Dim,&Eps,&Tau,&IterLim)<4 )
   {  cprintf("\n\r      Error reading parameters of the algorithm! \n");
      sound(120);  delay(800);  nosound();  goto final;
   }

if( Dim<1 )
   {  cprintf("\n\r      Incorrect dimension of the system! \n");
      sound(120);  delay(800);  nosound();  goto final;
   }

if( Eps<0 )
   {  cprintf("\n\r     Incorrect defect of the solution! \n");
      sound(120);  delay(800);  nosound();  goto final;
   }

DIM=2*Dim;
ip=calloc(DIM,sizeof(int));
G=calloc(DIM*DIM,sizeof(double));
x=calloc(DIM,sizeof(double));
xx=calloc(DIM,sizeof(double));
H=calloc(Dim*DIM,sizeof(double));
Z=calloc(Dim*DIM,sizeof(char));
d=calloc(DIM,sizeof(double));

if( ip==NULL||G==NULL||x==NULL||xx==NULL||H==NULL||Z==NULL||d==NULL )
   {  cprintf("\n\r    Not enough room to allocate the data! \n");
      sound(120);  delay(1000);  nosound();  goto final;
   }

for( i=0; i<Dim; i++ )
for( j=0; j<Dim; j++ )
    {  if( fscanf(fi,"%lf %lf",&u0,&u1)<2 )
	  {  cprintf("\n\r    Error reading the interval matrix! \n");
	     sound(120);  delay(800);  nosound();  goto final;
	  }
       *(H+i*DIM+2*j)=u0; *(H+i*DIM+2*j+1)=u1;
       p=0.5*(u0+u1);
       if( p>=0. )
	  { *(G+i*DIM+j)=p; *(G+(i+Dim)*DIM+j+Dim)=p;
	    *(G+(i+Dim)*DIM+j)=0; *(G+i*DIM+j+Dim)=0;
	  }
       else
	  { *(G+i*DIM+j)=0; *(G+(i+Dim)*DIM+j+Dim)=0;
	    *(G+(i+Dim)*DIM+j)=p; *(G+i*DIM+j+Dim)=p;
	  }
    }

for( i=0; i<Dim; i++ )
    {  if( fscanf(fi,"%lf %lf",&v0,&v1)<2 )
	  {  cprintf("\n\r    Error reading the right-hand side vector! \n");
	     sound(120);  delay(800);  nosound();  goto final;
	  }
       *(xx+i)=v0;  *(xx+i+Dim)=v1;
       *(d+2*i)=v0;  *(d+2*i+1)=v1;
    }

fclose(fi);

       /*   для контроля выводим на печать данные задачи      */

textcolor(LIGHTGRAY);
cprintf("\r\n Dimension of the system = %u",Dim);
cprintf("\r\n Interval matrix of the system: \n\n\r");

for( i=0; i<Dim; i++ )
 {  for( j=0; j<Dim; j++ )
	cprintf(" [ %- .3f, %- .3f ]",*(H+i*DIM+2*j),*(H+i*DIM+2*j+1));
    cprintf("\n\r");
 }

cprintf("\n Interval right-hand side vector: \n\n\r");

for( i=0; i<Dim; i++ )
    cprintf(" [ %- .3f, %- .3f ]",*(d+2*i),*(d+2*i+1));

cprintf("\n\n");  textcolor(RED+BLINK);
cprintf("\r%60s"," Wait!  The computer is working ... ");
textcolor(LIGHTGRAY);

   /*   solving the immersed "middle" system,
			  we find a starting approximation    */

     ier=0;     *(ip+DIM-1)=1;
for( k=0; k<DIM-1; k++ )
  { m=k; p=fabs(*(G+k*DIM+k));
    for( i=k+1; i<DIM; i++ )
      { q=fabs(*(G+i*DIM+k)); if( p<q ) { m=i; p=q; } }
      *(ip+k)=m; p=*(G+m*DIM+k);
      if( m!=k ) { *(ip+DIM-1)=-*(ip+DIM-1);
		   *(G+m*DIM+k)=*(G+k*DIM+k);
		   *(G+k*DIM+k)=p;
		 }
      if( p==0 ) { ier=k; *(ip+DIM-1)=0; break; }
      p=1./p;
      for( i=k+1; i<DIM; i++ ) *(G+i*DIM+k)*=-p;
      for( j=k+1; j<DIM; j++ )
	  { p=*(G+m*DIM+j); *(G+m*DIM+j)=*(G+k*DIM+j); *(G+k*DIM+j)=p;
	    if( p!=0 ) for( i=k+1; i<DIM; i++ ) *(G+i*DIM+j)+=*(G+i*DIM+k)*p;
	  }
  }

if(ier!=0 || *(G+DIM*DIM-1)==0)
    { delline(); sound(120); delay(1000); nosound(); textcolor(LIGHTGRAY);
      cprintf("\n\r   Unfortunately, the matrix of the system is i-singular!");
      cprintf("\n\r   Try to change its elements a little.\n");
      goto final;
    }

for( k=0; k<DIM-1; k++ )
 {  m=*(ip+k); q=*(xx+m); *(xx+m)=*(xx+k); *(xx+k)=q;
    for( i=k+1; i<DIM; i++ ) *(xx+i)+=*(G+i*DIM+k)*q;
 }

for( j=0; j<DIM-1; j++ )
  {   k=DIM-j-1;  *(xx+k)/=*(G+k*DIM+k);  q=-*(xx+k);
    for( i=0; i<k; i++ ) *(xx+i)+=*(G+i*DIM+k)*q;
  }
	     *xx/=*G; 


   /*      forming the splitting of the interval matrix     */

for( i=0; i<Dim; i++ )
for( j=0; j<Dim; j++ )
  { u0=*(H+DIM*i+2*j); u1=*(H+DIM*i+2*j+1); 
    p=0; *(Z+DIM*i+2*j)=(u0<u1)? 0:1;	   
     if( u0<=0 && u1<=0 ) p=(u0<u1)? u0:u1;       
     if( u0>=0 && u1>=0 ) p=(u0>u1)? u0:u1;
    *(H+DIM*i+2*j)=u0-p; *(H+DIM*i+2*j+1)=u1-p;
    *(Z+DIM*i+2*j+1)=(p==0)? 0:1;     
    if( p>=0. )
       { *(G+i*DIM+j)=p; *(G+(i+Dim)*DIM+j+Dim)=p;
	 *(G+(i+Dim)*DIM+j)=0; *(G+i*DIM+j+Dim)=0;
       }
    else
       { *(G+i*DIM+j)=0; *(G+(i+Dim)*DIM+j+Dim)=0;
	 *(G+(i+Dim)*DIM+j)=p; *(G+i*DIM+j+Dim)=p;
       }
  }

      /*    finding triangular decomposition of the splitted matrix G    */
	      
     ier=0;     *(ip+DIM-1)=1;
for( k=0; k<DIM-1; k++ )
  { m=k; p=fabs(*(G+k*DIM+k));
    for( i=k+1; i<DIM; i++ )
      { q=fabs(*(G+i*DIM+k)); if( p<q ) { m=i; p=q; } }
      *(ip+k)=m; p=*(G+m*DIM+k);
      if( m!=k ) { *(ip+DIM-1)=-*(ip+DIM-1);
		   *(G+m*DIM+k)=*(G+k*DIM+k);
		   *(G+k*DIM+k)=p;
		 }
      if( p==0 ) { ier=k; *(ip+DIM-1)=0; break; }
      p=1./p;
      for( i=k+1; i<DIM; i++ ) *(G+i*DIM+k)*=-p;
      for( j=k+1; j<DIM; j++ )
	  { p=*(G+m*DIM+j); *(G+m*DIM+j)=*(G+k*DIM+j); *(G+k*DIM+j)=p;
	    if( p!=0 ) for( i=k+1; i<DIM; i++ ) *(G+i*DIM+j)+=*(G+i*DIM+k)*p;
	  }
  }

if(ier!=0 || *(G+DIM*DIM-1)==0)
    { delline(); sound(120); delay(1000); nosound(); textcolor(LIGHTGRAY);
      cprintf("\n\r   Unfortunately, the matrix of the system is i-singular!");
      cprintf("\n\r   Try to change its elements a little.\n");
      goto final;
    }


   /*    launching iterations   */

             ni=0;

do {     ni++;   for( i=0; i<DIM; i++ ) *(x+i)=*(xx+i); 

      for( i=0; i<Dim; i++ )
       {    s0=*(d+2*i);   s1=*(d+2*i+1);
	  for( j=0; j<Dim; j++ )
	   { g0=*(H+i*DIM+2*j); g1=*(H+i*DIM+2*j+1); h0=*(x+j); h1=*(x+j+Dim);

	     if( *(Z+DIM*i+2*j+1) ) { p=h0; h0=h1; h1=p; }  
	     if( h0>=0 && h1>=0 )  m=1;
	     if( h0<=0 && h1>=0 )  m=2;  
	     if( h0<=0 && h1<=0 )  m=3;
	     if( h0>=0 && h1<=0 )  m=4;	

	     switch ( *(Z+DIM*i+2*j)*4+m ) 
		{ case 1: t0=g0*h1; t1=g1*h1; break;
		  case 2: p=g0*h1; q=g1*h0; t0=(p<q)?p:q; 
			  p=g0*h0; q=g1*h1; t1=(p>q)?p:q; break;
		  case 3: t0=g1*h0; t1=g0*h0; break;
		  case 4: t0=0; t1=0; break;
		  case 5: t0=g0*h0; t1=g1*h0; break;
		  case 6: t0=0; t1=0; break;
		  case 7: t0=g1*h1; t1=g0*h1; break;
		  case 8: p=g0*h0; q=g1*h1; t0=(p>q)?p:q; 
			  p=g0*h1; q=g1*h0; t1=(p<q)?p:q; break; 
		}
		   s0-=t0;       s1-=t1;
	   }
	      *(xx+i)=s0;     *(xx+i+Dim)=s1;  
	}

	/*      Solving a linear system with the matrix G.
		     Its triangular decomposition has been already found.   */
	    
      for( k=0; k<DIM-1; k++ )
	{  m=*(ip+k); q=*(xx+m); *(xx+m)=*(xx+k); *(xx+k)=q;
	  for( i=k+1; i<DIM; i++ ) *(xx+i)+=*(G+i*DIM+k)*q;
	}

      for( j=0; j<DIM-1; j++ )
	{  k=DIM-j-1;  *(xx+k)/=*(G+k*DIM+k);  q=-*(xx+k);
	  for( i=0; i<k; i++ ) *(xx+i)+=*(G+i*DIM+k)*q;
	}
		   *xx/=*G; 

	/*    computing relative improvement of the current
					approximation to the solution    */

     p=q=0; 
     for( i=0; i<DIM; i++ ) { q+=fabs(*(xx+i)); p+=fabs(*(xx+i)-*(x+i)); }
     if( q==0 ) q=1;

   } while ( p/q>Eps && ni<IterLim );

delline();
if( ni<IterLim )
  { sound(600); delay(150); sound(400); delay(150); sound(600); delay(300);
    nosound(); delay(50); sound(400); delay(300); nosound();
  }
else
  { sound(100); delay(1500); nosound(); textcolor(LIGHTGRAY);
    cprintf("\r       Doesn't the method diverge?! \n");
    cprintf("\r       The number of iterations alloted has been exhausted!\n");
  }
cprintf("\r Algebraic solution of the interval system:\n");

for( i=0; i<Dim; i++ )
    {  u0=*(xx+i);  u1=*(xx+i+Dim);
       cprintf("\r\n%15d%7s",i+1,".     ["); textcolor(WHITE);
       cprintf(" %16.11f",u0); textcolor(LIGHTGRAY);
       cprintf(" ,"); textcolor(WHITE); cprintf("%16.11f",u1);
       textcolor(LIGHTGRAY); cprintf(" ]");
       if( u0<=u1 )  cprintf("     ->");
       else  cprintf("     <-");
    }

cprintf("\n\n\r Number of iterations = %d",ni);
cprintf("\n\r 1-norm of the defect of the approximate solution = %.3E\n",p);

final: textcolor(RED);
cprintf("\n\r %48s","Press a key to exit ");
getch(); delline(); textcolor(LIGHTGRAY); cprintf("\r %40s","Bye-bye!");
delay(500); delline();

}