/*     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 <string.h>

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

    printf("\r\n\n     *** ");
    printf("  Computing  formal  solutions  to  interval  linear  systems  ");
    printf(" *** \r\n      ** ");
    printf("  by iterative method with real splitting, (c) S. Shary, 1995  ");
    printf(" ** \n");

    if (argc == 1){
        printf("\n\rProgram 'resplit'  computes  formal  solutions  to  interval");
        printf(" linear systems  in\n");
        printf("\r Kaucher arithmetic  and implements stationary one-step");
        printf(" iterative method  based\n");
        printf("\r on  modified  real splitting  of the matrix  of the system,");
        printf(" which  is supposed\n");
        printf("\r to be square. \n\n");
        printf("\rUsage:");
        printf("   resplit  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("an auxiliary real parameter (may be arbitrary),\n");
        printf("\r %15s"," * ");
        printf("restriction to the number of iterations (a positive integer). \n\n");
        printf("\rNext the interval matrix (along the rows) and");
        printf(" right-hand side vector are written out in the DataFile.\n\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);
        goto final;
    }

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

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

    if (Eps < 0){
        printf("\n\r     Incorrect defect of the solution! \n");
        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)){
        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;
            }
            *(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){
            printf("\n\r    Error reading the right-hand side vector! \n");
            goto final;
        }
        *(xx+i)=v0;
        *(xx+i+Dim)=v1;
        *(d+2*i)=v0;
        *(d+2*i+1)=v1;
    }

    fclose(fi);

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

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

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

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

    printf("\n\n");
    printf("\r%60s"," Wait!  The computer is working ... ");
    
    /*   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)){
        printf("\n\r   Unfortunately, the matrix of the system is i-singular!");
        printf("\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)){
        printf("\n\r   Unfortunately, the matrix of the system is i-singular!");
        printf("\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));

    if(ni >= IterLim){
        printf("\r       Doesn't the method diverge?! \n");
        printf("\r       The number of iterations alloted has been exhausted!\n\n");
    }
    printf("\r Algebraic solution of the interval system:\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",p);
    final: ; 
}