/***************************************************************************
*          Solve a symmetric linear system using function SYMSOL           *
* ------------------------------------------------------------------------ *
* SAMPLE RUN:                                                              *
*  ( Solve symmetric linear system:                                        *
*      X1 -2X2 +3X3 -  4X4 +     5X5 -     6X6 = -3                        *
*    -2X1 +3X2 -5X3 +  7X4 -     4X5 +     8X6 =  7                        *
*     3X1 -5X2 +5X3 +  6X4 -     2X5 +      X6 =  8                        *
*    -4X1 +7X2 +6X3 -  9X4 +    20X5 +   0.5X6 = 20.5                      *
*     5X1 -4X2 -2X3 + 20X4 +    11X5 +0.3333X6 = 30.3333                   *
*    -6X1 +8X2 + X3 +0.5X4 +0.3333X5 +    13X6 = 16.8333 )                 *
*                                                                          *
*  SYMMETRIC SYSTEM TO SOLVE:                                              *
*  1.000000 -2.000000  3.000000 -4.000000  5.000000  -6.000000   -3.000000 *
* -2.000000  3.000000 -5.000000  7.000000 -4.000000   8.000000    7.000000 *
*  3.000000 -5.000000  5.000000  6.000000 -2.000000   1.000000    8.000000 *
* -4.000000  7.000000  6.000000 -9.000000 20.000000   0.500000   20.500000 *
*  5.000000 -4.000000 -2.000000 20.000000 11.000000   0.333300   30.333300 *
* -6.000000  8.000000  1.000000  0.500000  0.333300  13.000000   16.833300 *
*                                                                          *
*  SOLUTIONS:                                                              *
*  1.00000000                                                              *
*  1.00000000                                                              *
*  1.00000000                                                              *
*  1.00000000                                                              *
*  1.00000000                                                              *
*  1.00000000                                                              *
*                                                                          *
* ------------------------------------------------------------------------ *
* Reference:  From Numath Library By Tuan Dang Trong in Fortran 77         *
*             [BIBLI 18].                                                  *
*                                                                          *
*                                  C++ Release 1.0 By J-P Moreau, Paris    *
*                                            (www.jpmoreau.fr)             *
****************************************************************************
To link with basis_r.cpp and vmblock.cpp. 
----------------------------------------- */
#include <basis.h>
#include <vmblock.h>


      REAL **A;
      REAL *B,*X,*D;
      int I,J,N;
      void *vmblock = NULL;


    void SYMSOL(int N, REAL **A, REAL *B, REAL *X, REAL *D)  {
/*---------------------------------------------------------------------
!     SOLVE A SYMMETRIC LINEAR SYSTEM A.X = B
!     A IS SUPPOSED WELL BALANCED AND CAN BE NOT POSITIVE DEFINITE,
!     THE COEFFICIENTS OF A ARE STORED IN THE UPPER TRIANGULAR HALF
!     OF THE MATRIX AND ARE LOST DURING THE PROCESS.
!     THE USED METHOD IS BASED UPON THE DECOMPOSITION OF THE INITIAL 
!     MATRIX INTO A PRODUCT OF 3 MATRICES: A = (U TRANSPOSE).D.U
!     INPUTS:
!     N       NUMBER OF LINES OF THE LINEAR SYSTEM
!     A       TABLE OF DIMENSION N X N STORING THE COEFFICIENTS
!             OF THE SYSTEM
!     B       VECTOR OF DIMENSION N (SECOND MEMBER)
!     OUTPUT:
!     X       SYSTEM SOLUTION
!     WORKING ZONE:
!     D       DIAGONAL D OF DECOMPOSITION
!---------------------------------------------------------------------*/
      // Labels: e2, e4, e6, e7, e8, e10
      REAL S;
      int  I,K,L,M;
      for (I=1; I<=N; I++) {
        L = I-1;
        M = I+1;
        S = A[I][I];
        if (I==1) goto e2;
        for (K=1; K<=L; K++)
          S = S - A[K][I]*D[K]*A[K][I];
e2:     D[I] = S;
        if (M > N) goto e6;
        for (J=M; J<=N; J++) {
          S = A[I][J];
          if (I==1) goto e4;
          for (K=1; K<=L; K++)
            S = S - A[K][I]*D[K]*A[K][J];
e4:       A[I][J] = S/D[I];
        }
      }
e6:   for (I=1; I<=N; I++) {
        S = B[I];
        L = I-1;
        if (I==1) goto e8;
        for (K=1; K<=L; K++)
          S = S - A[K][I]*X[K];
e8:     X[I] = S;
      }
      for (L=1; L<=N; L++) {
        I = N-L+1;
        S = X[I]/D[I];
        M = I+1;
        if (M>N) goto e10;
        for (K=M; K<=N; K++)
          S = S - A[I][K]*X[K];
e10:    X[I] = S;
      }
    } //SYMSOL()

void main()  {

  N = 6;    // size of linear system

  // allocate memory for matrix A and vectors B, X, D (index 0 not used)
  vmblock = vminit();
  A   = (REAL **)vmalloc(vmblock, MATRIX, N+1, N+1);
  X   = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
  B   = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
  D   = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);

  if (! vmcomplete(vmblock))
    LogError ("No Memory", 0, __FILE__, __LINE__);

  //define upper half part of matrix A
  A[1][1]= 1.0; A[1][2]=-2.0; A[1][3]=3.0; A[1][4]=-4.0; A[1][5]=5.0; A[1][6]=-6.0;
  A[2][2]= 3.0; A[2][3]=-5.0; A[2][4]=7.0; A[2][5]=-4.0; A[2][6]=8.0;                
  A[3][3]= 5.0; A[3][4]= 6.0; A[3][5]=-2.0; A[3][6]=1.0;
  A[4][4]=-9.0; A[4][5]=20.0; A[4][6]=0.50;
  A[5][5]=11.0; A[5][6]= 0.3333;
  A[6][6]=13.0;

  //define 2nd member
  B[1]=-3.0;
  B[2]=7.0;
  B[3]=8.0;
  B[4]=20.5;
  B[5]=30.3333;
  B[6]=16.8333;

  // optional section to print full data
  for (I=2; I<=N; I++)
    for (J=1; J<I; J++)
      A[I][J] = A[J][I];

  printf("\n SYMMETRIC SYSTEM TO SOLVE:\n");
  for (I=1; I<=N; I++) {
    for (J=1; J<=N; J++) printf("%10.6f", A[I][J]);
    printf("  %10.6f\n", B[I]);
  }
  // end optional section

  // call solver for symmetric linear systems
  SYMSOL(N,A,B,X,D);

  //print solutions
  printf(" SOLUTIONS:\n");
  for (I=1; I<=N; I++) printf("%12.8f\n", X[I]);

  // free memory
  vmfree(vmblock);
  printf("\n");

}

// end of file tsymsol.cpp