/*       Computing formal-algebraic solution to interval linear systems 
	 of the form x=Cx+b in Kaucher arithmetic by subddifferential 
	      Newton method with a special starting approximation, 
	   (c) Sergey Shary, 1996                                          */

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

     int   DIM, *ip;
  double   *C, *d, *F, *x, *xx;


char solver()

   /*    solving a point linear system of the double 
			      dimension by Gaussian elimination    */
{     int   ier, i, j, k, m;
   double   p, q;

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

  if(ier!=0 || *(F+DIM*DIM-1)==0)
      { delline(); sound(120); delay(1000); nosound(); textcolor(LIGHTGRAY);
	cprintf("\n\r   Sorry, the matrix of the equation is i-singular."); 
	cprintf("\n\r   Try to change it a little. \n"); 
	return 1;
      }

  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)+=*(F+i*DIM+k)*q;
   }

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


main(int argc, char *argv[])

{   char   DataFile[64];
     int   i, j, l, m, ni, Dim;
  double   g0, g1, h0, h1, p, q, r, s0, s1, t0, t1, u0, u1, v0, v1, Eps, Tau;
    FILE   *fi, *fopen();
    long   IterLim;

textmode(C80); clrscr(); textcolor(RED); cprintf("\r\n  ***   ");
textcolor(LIGHTRED);
cprintf("Computing formal-algebraic solutions of interval linear equations ");
textcolor(RED); cprintf("  *** \r\n   ** "); textcolor(LIGHTRED);
cprintf("    of the form  x = Cx + b  by  subdifferential  Newton  method ");
textcolor(RED); cprintf("    ** \r\n    * "); textcolor(LIGHTRED);
cprintf("                     (c) Sergey P. Shary, 1995                ");
textcolor(RED); cprintf("       * \n\n"); delay(1200); textcolor(LIGHTGRAY);

if(argc == 1)
{ cprintf("\n\rThe ALOUTER0 program  implements  subdifferential Newton ");
  cprintf("method and is intended");
  cprintf(" for  the computation  of  formal  solutions  to  interval ");
  cprintf("linear equations of the form  x = Cx + b  in");
  cprintf(" Kaucher complete arithmetic.  The interval matrix C ");
  cprintf("is supposed to be square.\n\n");
  cprintf("\rUsage:"); textcolor(WHITE); 
  cprintf("   ALOUTER0  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("damping factor (a real within (0,1], 1. is recommended),\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 C (along the rows) and");
  cprintf(" vector b  are written out in the DataFile.\n"); 
  goto final;
}
else
  strcpy(DataFile,argv[1]); 

  /*  reading input data of the problem,
		  forming the "middle" point matrix F 
				  of the interval system under solution    */

                          fi=fopen(DataFile,"r");
if( fi==NULL )
{  printf("\n\r    Cannot open the data file `%1s' \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 the parameters of the algorithm! \n");  
      sound(120);  delay(800);  nosound();  goto final;
   }

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

if( Tau<=0 || Tau>1 )
   {  cprintf("\n\r      Incorrect damping factor! \n"); 
      sound(120);  delay(800);  nosound();  goto final;
   }

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

DIM=2*Dim; 
ip=calloc(DIM,sizeof(int));
F=calloc(DIM*DIM,sizeof(double));
x=calloc(DIM,sizeof(double));
xx=calloc(DIM,sizeof(double));
C=calloc(Dim*DIM,sizeof(double));
d=calloc(DIM,sizeof(double));

if( ip==NULL || F==NULL || x==NULL || xx==NULL || C==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;
	  }
       if( u0>u1 )
	  {  cprintf("\n\rThe entry (%d,%d) of the interval matrix",i+1,j+1);
	     cprintf("\n is improper! \n");
	     sound(120);  delay(800);  nosound();  goto final;
	  }  
       *(C+i*DIM+2*j)=u0; *(C+i*DIM+2*j+1)=u1;
       p=-0.5*(u0+u1); if( i==j ) p+=0.5;
       if( p>=0. )
	  { *(F+i*DIM+j)=p; *(F+(i+Dim)*DIM+j+Dim)=p;
	    *(F+(i+Dim)*DIM+j)=0; *(F+i*DIM+j+Dim)=0; 
	  }
       else
	  { *(F+i*DIM+j)=0; *(F+(i+Dim)*DIM+j+Dim)=0;
	    *(F+(i+Dim)*DIM+j)=p; *(F+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;
	  }
       if( v0>v1 ) 
	  {  cprintf("\n\rThe %d-th entry of the right-side vector",i+1);
	     cprintf("\n is improper!\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);

       /*     displaying the input data of the problem      */

textcolor(LIGHTGRAY);
cprintf("\r\n Dimension of the equation = %u",Dim);
cprintf("\r\n The interval matrix C  of the equation: \n\n\r");

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

cprintf("\n The interval vector b: \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%57s"," Wait! The computer is working ... ");
textcolor(LIGHTGRAY);

    /*   solving the embedded "middle" point system,
		   we find the starting iteration of the algorithm   */

if( solver() ) goto final;

   /*    launching the iteration 
			 of the subdifferential Newton method    */
             ni=0;

do {  ni++; r=0;  for( i=0; i<DIM; i++ )
		      {  *(x+i)=*(xx+i);  
			 for( j=0; j<DIM; j++ ) *(F+i*DIM+j)= (i!=j)? 0:1;
		      }

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

             if ( g0*g1>0 ) l= (g0>0)? 0:2; else l= (g0<=g1)? 1:3;
             if ( h0*h1>0 ) m= (h0>0)? 1:3; else m= (h0<=h1)? 2:4;

             switch ( 4*l+m ) {
             case  1: t0=g0*h0; t1=g1*h1;
		      *(F+i*DIM+j)-=g0; *(F+(i+Dim)*DIM+j+Dim)-=g1;
		      break;
	     case  2: t0=g1*h0; t1=g1*h1;
		      *(F+i*DIM+j)-=g1; *(F+(i+Dim)*DIM+j+Dim)-=g1;
		      break;
	     case  3: t0=g1*h0; t1=g0*h1;
		      *(F+i*DIM+j)-=g1; *(F+(i+Dim)*DIM+j+Dim)-=g0;
		      break;
	     case  4: t0=g0*h0; t1=g0*h1;
		      *(F+i*DIM+j)-=g0; *(F+(i+Dim)*DIM+j+Dim)-=g0;
		      break;
	     case  5: t0=g0*h1; t1=g1*h1;
		      *(F+i*DIM+j+Dim)-=g0; *(F+(i+Dim)*DIM+j+Dim)-=g1;
		      break;
	     case  6: u0=g0*h1; v0=g1*h0;  u1=g0*h0; v1=g1*h1;
		      if ( u0<v0 ) { t0=u0; *(F+i*DIM+j+Dim)-=g0; }
		      else { t0=v0; *(F+i*DIM+j)-=g1; }
		      if ( u1>v1 ) { t1=u1; *(F+(i+Dim)*DIM+j)-=g0; }
		      else { t1=v1; *(F+(i+Dim)*DIM+j+Dim)-=g1; }
		      break;
	     case  7: t0=g1*h0; t1=g0*h0;
		      *(F+i*DIM+j)-=g1; *(F+(i+Dim)*DIM+j)-=g0;
		      break;
	     case  8: t0=0; t1=0;
		      break;
	     case  9: t0=g0*h1; t1=g1*h0;
		      *(F+i*DIM+j+Dim)-=g0; *(F+(i+Dim)*DIM+j)-=g1;
		      break;
	     case 10: t0=g0*h1; t1=g0*h0;
		      *(F+i*DIM+j+Dim)-=g0; *(F+(i+Dim)*DIM+j)-=g0;
		      break;
	     case 11: t0=g1*h1; t1=g0*h0;
		      *(F+i*DIM+j+Dim)-=g1; *(F+(i+Dim)*DIM+j)-=g0;
		      break;
	     case 12: t0=g1*h1; t1=g1*h0;
		      *(F+i*DIM+j+Dim)-=g1; *(F+(i+Dim)*DIM+j)-=g1;
		      break;
	     case 13: t0=g0*h0; t1=g1*h0;
		      *(F+i*DIM+j)-=g0; *(F+(i+Dim)*DIM+j)-=g1;
		      break;
	     case 14: t0=0; t1=0;
		      break;
	     case 15: t0=g1*h1; t1=g0*h1;
		      *(F+i*DIM+j+Dim)-=g1; *(F+(i+Dim)*DIM+j+Dim)-=g0;
		      break;
	     case 16: u0=g0*h0; v0=g1*h1;  u1=g0*h1; v1=g1*h0;
		      if (u0>v0) { t0=u0; *(F+i*DIM+j)-=g0; }
		      else { t0=v0; *(F+i*DIM+j+Dim)-=g1; }
		      if (u1<v1) { t1=u1; *(F+(i+Dim)*DIM+j+Dim)-=g0; }
		      else { t1=v1; *(F+(i+Dim)*DIM+j)-=g1; }
		      break;  }
		   s0-=t0;       s1-=t1;
	   }
	     *(xx+i)=t0=s0-*(d+2*i);   *(xx+i+Dim)=t1=s1-*(d+2*i+1);
	     t0=fabs(t0); t1=fabs(t1);     r+= (t0>t1)? t0:t1;
       }

		/*      solving the point linear system 
				       with the subgradient matrix F     */

     if( solver() ) goto final;

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

     if( q==0 ) q=1;

   } while ( r/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 is exhausted!\n");
  }
cprintf("\r Formal-algebraic solution of the equation:\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",r);

final: textcolor(RED);
cprintf("\n\r %50s","Press any key to quit "); 
getch(); delline(); textcolor(LIGHTGRAY); cprintf("\r %42s","Bye-bye!"); 
delay(500); delline();  

}