/*         Computing formal-algebraic solution  to interval linear systems 
	 of the form Gx=(G-A)x+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   *G, *H, *b, *F, *x, *xx;

	  /*  Global variables:

	      DIM - double dimension of the original problem, which is 
		    actually the dimension of the auxiliary point equation
		    to be solved in the embedded Euclidean space. 
	      ip  - auxiliary integer array that is necessary for Gaussian
		    Gaussian elimination in the function "solver". It keeps 
		    information about column permutations.  			      
	      G   - pointer to the storage area where the point diagonal 
		    matrix splitted from A  is allocated
	      H   - pointer to the storage area where the interval matrix
		    (G-A) is allocated
	      b   - pointer to the storage area where the right-hand side
		    interval vector is allocated        		   
	      F	  - pointer to the storage area where the auxiliary matrix
		    of the double size is allocated. During the iteration
		    of the subdifferential Newton method, F is equal to the
		    the subgradient of the iterated function (in the embedded
		    double dimension space)                                 */


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;

	 /*      displaying the title of the program    */
	   
textmode(C80); clrscr(); textcolor(RED); cprintf("\r\n   ***   ");
textcolor(LIGHTRED);
cprintf("Solving the `outer problem' for  interval  linear  system  ");
textcolor(RED); cprintf("  *** \r\n    ** "); textcolor(LIGHTRED);
cprintf("  of the form   Ax = 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);

	 /*      displaying information about itself if need be,
			 reading the name of the datafile otherwise      */

if(argc == 1)
{ cprintf("\n\rThe ALOUTER2 program is intended for the solution of the ");
  cprintf("`outer problem' for the interval linear systems of the form ");
  cprintf(" Ax = b  with the square interval matrix A. "); 
  cprintf("ALOUTER2 relyes on the 'formal' approach to the `outer problem'");
  cprintf(" and implements the subdifferential Newton method");
  cprintf(" applied to the fixed point equation, which is obtained by");
  cprintf(" splitting the deviation of the main diagonal of A.\n\n");
  cprintf("\rUsage:"); textcolor(WHITE); 
  cprintf("   ALOUTER2  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 A (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        */

                          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;
   }

	     /*   allotting the storage areas for input data  */ 

DIM=2*Dim; 
ip=calloc(DIM,sizeof(int));
F=calloc(DIM*DIM,sizeof(double));
x=calloc(DIM,sizeof(double));
xx=calloc(DIM,sizeof(double));
G=calloc(Dim,sizeof(double));
H=calloc(Dim*DIM,sizeof(double));
b=calloc(DIM,sizeof(double));

if( ip==NULL||F==NULL||x==NULL||xx==NULL||G==NULL||H==NULL||b==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(" is improper! \n");
	     sound(120);  delay(800);  nosound();  goto final;
	  }  
       *(H+i*Dim+j)=u0; *(H+(i+Dim)*Dim+j)=u1;
    }

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(" is improper!\n");
	     sound(120);  delay(800);  nosound();  goto final;
	  }
       *(b+i)=v0;  *(b+i+Dim)=v1;
       *(xx+i)=0;  *(xx+i+Dim)=0;
    }

fclose(fi);

	/*    displaying the input data of the problem, 
				     forming the additives G and H      */

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

for( i=0; i<Dim; i++ )
 {  for( j=0; j<Dim; j++ )
	{  cprintf(" [ %- .3f, %- .3f ]",p=*(H+i*Dim+j),q=*(H+(i+Dim)*Dim+j));
	   if( i==j ) 
	      {  r= (fabs(p) > fabs(q))? p:q;
		 *(G+i)=r; *(H+i*Dim+i)=r-q; *(H+(i+Dim)*Dim+i)=r-p; 
	      }
	    else
	      { *(H+i*Dim+j)=-q; 
		*(H+(i+Dim)*Dim+j)=-p;
	      }

	}
    cprintf("\n\r");
 }

cprintf("\n The interval right-hand side vector b: \n\n\r");

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

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


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

do {    ni++;      r=0;     

      for( i=0; i<DIM; i++ )  
       {  *(x+i)=*(xx+i);  
	 for( j=0; j<DIM; j++ )  *(F+i*DIM+j)=0;
       }

      for( i=0; i<Dim; i++ )
       {     
		      p=*(G+i); 

	     if( p>=0. )
		{ *(F+i*DIM+i)+=p; *(F+(i+Dim)*DIM+i+Dim)+=p;
		  s0=*(x+i)*p; s1=*(x+i+Dim)*p; 
		}
	     else
		{ *(F+(i+Dim)*DIM+i)+=p; *(F+i*DIM+i+Dim)+=p; 
		  s0=*(x+i+Dim)*p; s1=*(x+i)*p; 
		} 

	 for( j=0; j<Dim; j++ )
	   { 
		g0=*(H+i*Dim+j);   g1=*(H+(i+Dim)*Dim+j); 
		      h0=*(x+j);   h1=*(x+j+Dim);

	       /*   identifying the types of the intervals H_{ij} and x_j   */

             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;

                /* computing multiplication in Kaucher arithmetic,
                                        forming the subgradient matrix F    */
             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;
	   }

	     t0=fabs( *(xx+i)=s0-*(b+i) ); 
	     t1=fabs( *(xx+i+Dim)=s1-*(b+i+Dim) );
	     r+= (t0>t1)? t0:t1;
       }  

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

     if( solver() ) goto final;

		/*      computing the next approximation 
				     to the solution  and its 1-norm    */

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

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


		 /*      concluding music      */
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");
  }


		 /*   outputting the answer  
					 and related parameters   */   
cprintf("\r Formal-algebraic solution of the equation");
cprintf("  (dev diag A) x = (dev diag A - A) x + b :\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);

		 /*           farewell           */

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

}