/*         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 <string.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)){
        printf("\n\r   Sorry, the matrix of the equation is i-singular.");
        printf("\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;
}



void 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    */

    printf("\r\n      ***   ");
    printf("    Solving the `outer problem' for  interval  linear  system ");
    printf("      *** \r\n       **   ");
    printf("    of the form   Ax = b  by  subdifferential  Newton  method");
    printf("%9.3s\r\n","**");  printf("%9.3s","*"); 
    printf("                       (c) Sergey P. Shary, 1995           ");
    printf("%13.1s\n\n","*");

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

    if (argc == 1){
        printf("\n\rThe program 'alouter2' is intended for the solution of the `outer problem' for the interval\n");
        printf("linear systems of the form Ax = b  with the square interval matrix A. 'alouter' relies\n");
        printf("on the 'formal' approach to the `outer problem' and implements the subdifferential Newton\n");
        printf("method applied to the fixed point equation, which is obtained by splitting the deviation\n"); 
        printf("of the main diagonal of A.\n\n");
        printf("\rUsage:");
        printf("   alouter2  DataFile \n\n ");
        printf("\rIn the DataFile, the following parameters  must occupy  the");
        printf(" first four positions (in order): \n");
        printf("\r %15s"," * ");
        printf("dimension of the interval equation (a positive integer), \n");
        printf("\r %15s"," * ");
        printf("relative defect of the approximate solution (a positive real),\n");
        printf("\r %15s"," * ");
        printf("damping factor (a real within (0,1], 1. is recommended),\n");
        printf("\r %15s"," * ");
        printf("restriction to the number of iterations (a positive integer). \n\n");
        printf("\rNext the interval matrix A (along the rows) and");
        printf(" vector b  are written out in the DataFile.\n\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);
        goto final;
    }

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

    if (Dim < 1)  {
        printf("\n\r      Incorrect dimension of the equation! \n");
        goto final;
    }

    if ((Tau <= 0) || (Tau > 1)){
        printf("\n\r      Incorrect damping factor! \n");
        goto final;
    }

    if (Eps < 0){
        printf("\n\r      Incorrect relative defect! \n");
        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)){
        printf("\n\r    Not enough room to allocate the data! \n");
        goto final;
    }

    for (i = 0; i < Dim; i++){
        for (j = 0; j < Dim; j++){
            if (fscanf(fi,"%lf %lf",&u0,&u1) < 2){
                printf("\n\r    Error reading the interval matrix! \n");
                goto final;
            }
            if (u0 > u1){
                printf("\n\rThe entry (%d,%d) of the interval matrix",i+1,j+1);
                printf(" is improper! \n");
                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){
            printf("\n\r    Error reading the right-hand side vector! \n");
            goto final;
        }
        if (v0 > v1){
            printf("\n\rThe %d-th entry of the right-side vector",i+1);
            printf(" is improper!\n");
            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      */


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

    for (i = 0; i < Dim; i++){
        for (j = 0; j < Dim; j++){
            printf(" [ %- .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;
            }
        }
        printf("\n\r");
    }

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

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

    printf("\n\n");
    printf("\r%57s"," Wait! The computer is working ... ");


   /*    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));

    if( ni<IterLim ){
    }
    else{
        printf("\r       Doesn't the method diverge?! \n");
        printf("\r       The number of iterations alloted is exhausted!\n");
    }


         /*   outputting the answer
                     and related parameters   */
    printf("\r Formal-algebraic solution of the equation");
    printf("  (dev diag A) x = (dev diag A - A) x + b :\n");

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

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

         /*           farewell           */

    final: ;
}