/*****************************************************************
*              Test Gauss method for band matrix                 *
* -------------------------------------------------------------- *
* Reference:                                                     *
*                                                                *
* "Numerical Algorithms with C, By Gisela Engeln-Muellges        *
*  and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].           *
* -------------------------------------------------------------- *
* SAMPLE RUN:                                                    *
*                                                                *
* Input file: bandmat.dat                                        *
*                                                                *
* 5                                                              * 
* 1                                                              *
* 1                                                              *
*  0.0   1.0   -2.0                                              *
* 10.0   5.0   -7.5                                              *
*  2.22 -3.25   0.00456                                          *
*  2.0   4.0   -3.3                                              *
* -1.0   3.0    0.0                                              *
*                                                                *
* Output file: bandmat.lst                                       *
*                                                                *
* -------------------------------------------------------------- *
* Banded matrix                                                  *
* -------------------------------------------------------------- *
*                                                                *
* Dimension     : 5                                              *
* Subdiagonals  : 1                                              *
* Superdiagonals: 1                                              *
*                                                                *
* Input condensed band matrix:                                   *
*                                                                *
*   0.000000     1.000000    -2.000000                           *
*  10.000000     5.000000    -7.500000                           *
*   2.220000    -3.250000     0.004560                           *
*   2.000000     4.000000    -3.300000                           *
*  -1.000000     3.000000     0.000000                           *
*                                                                *
* Input uncondensed band matrix                                  *
*                                                                *
*   1.000000    -2.000000     0.000000     0.000000     0.000000 *
*  10.000000     5.000000    -7.500000     0.000000     0.000000 *
*   0.000000     2.220000    -3.250000     0.004560     0.000000 * 
*   0.000000     0.000000     2.000000     4.000000    -3.300000 *
*   0.000000     0.000000     0.000000    -1.000000     3.000000 *
*                                                                *
* -------------------------------------------------------------- *
* Band without pivot:                                            *
* -------------------------------------------------------------- *
*                                                                *
* Inverse (transposed)                                           *
*                                                                *
*  -0.005941    -0.502971    -0.343236     0.236714     0.078905 *
*   0.100594     0.050297     0.034324    -0.023671    -0.007890 *
*  -0.231916    -0.115958    -0.386526     0.266570     0.088857 *
*   0.000365     0.000182     0.000608     0.344408     0.114803 *
*   0.000401     0.000201     0.000669     0.378849     0.459616 *
*                                                                *
* Determinant = -5.627040e+002                                   *
*                                                                *
* -------------------------------------------------------------- *
* Band with pivot                                                *
* -------------------------------------------------------------- *
*                                                                *
* Inverse (transposed)                                           * 
*                                                                *
*  -0.005941    -0.502971    -0.343236     0.236714     0.078905 * 
*   0.100594     0.050297     0.034324    -0.023671    -0.007890 *
*  -0.231916    -0.115958    -0.386526     0.266570     0.088857 *
*   0.000365     0.000182     0.000608     0.344408     0.114803 *
*   0.000401     0.000201     0.000669     0.378849     0.459616 *
*                                                                *
* Determinant =  5.627040e+002                                   *
*                                                                *
* System Solution                                                * 
*                                                                *
*  -0.136497    -0.568249    -0.694162     1.202870     0.734290 *
*                                                                *
* -------------------------------------------------------------- *
* Uses files:                                                    *                                                        *
*              basis_r.cpp, vmblock.cpp, fband.cpp, fbando.cpp   *
*****************************************************************/
#include "basis.h"
#include "vmblock.h" 

//headers of called functions   
int bando(int, int, int, int, REAL **, REAL *); 
int band (int, int, int, int, REAL **, REAL *, int *, int *);

#define   SIZE  25                                                         


int main ()  {

  int ud, ld, n, i, k, rc, *perm, dsign, dim, mod;
  REAL **a, **packmat, **save, *b, determ;
  void *vmblock = NULL; 
   
  REAL solution[SIZE];
  REAL matinv[SIZE][SIZE];
  FILE *fp1, *fp2;

  fp1=fopen("bandmat5.dat","r");   
  fp2=fopen("bandmat5.lst","w");        
           
  WriteHead(fp2," Banded matrix");

  if (fscanf(fp1,"%d", &n) <= 0)
  {
    LogError ("Input stream", 0,  __FILE__, __LINE__);
    return 1;
  }

  if (n < 1)
  {
    LogError ("n must be > 0", 0,  __FILE__, __LINE__);
    return 1;
  }

  if (fscanf(fp1,"%d", &ld) <= 0)
  {
    LogError ("Input stream", 0,  __FILE__, __LINE__);
    return 1;
  }

  if (ld < 0)
  {
    LogError ("ld must be > 0", 0,  __FILE__, __LINE__);
    return 1;
  }

  if (fscanf(fp1,"%d", &ud) <= 0)
  {
    LogError ("Input stream", 0,  __FILE__, __LINE__);
    return 1;
  }

  if (ud < 0)
  {
    LogError ("ud must be > 0", 0,  __FILE__, __LINE__);
    return 1;
  }

  fprintf (fp2,"Dimension     : %d\n",      n );
  fprintf (fp2,"Subdiagonals  : %d\n",   ld);
  fprintf (fp2,"Superdiagonals: %d\n\n", ud);

  if (ld + ud >= n)
  {
    LogError ("ld + ud must be < n", 0, __FILE__, __LINE__);
    return (1);
  }

  dim = ld + ud + 1 + min (ld, ud);

  vmblock = vminit();
  packmat = (REAL **)vmalloc(vmblock, MATRIX,  n, dim);   
  a       = (REAL **)vmalloc(vmblock, MATRIX,  n, n);  
  save    = (REAL **)vmalloc(vmblock, MATRIX,  n, dim);
  b       = (REAL *) vmalloc(vmblock, VEKTOR,  n, 0);
  perm    = (int *)  vmalloc(vmblock, VVEKTOR, n, sizeof(*perm));

  if (! vmcomplete(vmblock))
  {
    LogError ("No Memory", 0, __FILE__, __LINE__);
    return 1;
  }
    
  if (ReadMat (fp1, n, ld + ud + 1, packmat))
  {
    LogError ("Input Stream", 0, __FILE__, __LINE__);
    return 1;
  }                
  fclose(fp1);

  CopyMat (n, ld + ud + 1, packmat, save);

  fprintf(fp2," Input condensed band matrix:\n");

  WriteMat(fp2, n, ld + ud + 1, packmat);  

  for (i=0; i<n; i++)
    for (k=0; k<n; k++)      
      if ((ld+k-i>-1) && (ld+k-i<ld+ud+1))
        a[i][k]=packmat[i][ld+k-i];  
      else
        a[i][k]=ZERO;  
      
  fprintf(fp2," Input uncondensed band matrix:\n"); 
  
  WriteMat (fp2, n, n, a);    

  WriteHead(fp2," Band without pivot:");
  fprintf (fp2," Inverse (transposed)\n\n");

  mod = 0;
  for (i = 0; i < n; i++)
  {
    SetVec (n, b, ZERO);
    b[i] = ONE;

    rc = bando (mod, n, ld, ud, packmat, b);
    if (rc)
    {
      LogError ("bando", rc, __FILE__, __LINE__);
      return 1;
    }
    WriteVec (fp2, n, b);
    mod = 2;
  }

  determ = ONE;
  for (i = 0; i < n; i++) determ *= packmat[i][ld];

  fprintf (fp2,"\n Determinant = ");
  fprintf (fp2,FORMAT_LE, determ);  
  fprintf(fp2,"\n");   

  WriteHead(fp2," Band with pivot");

  fprintf(fp2," Inverse (transposed)\n\n");

  mod = 0;
  for (i = 0; i < n; i++)
  {
    SetVec (n, b, ZERO);
    b[i] = ONE;

    rc = band (mod, n, ld, ud, save, b, perm, &dsign);
    if (rc)
    {
      LogError ("band", rc, __FILE__, __LINE__);
      return 1;
    }       
    WriteVec(fp2, n, b); 
    
    if (n<=SIZE) {                     
      solution[i]=ZERO;
      for (k=0; k<n; k++) matinv[i][k]=b[k];  
    }  
    mod = 2;
  }

  determ = (REAL) dsign;
  for (i = 0; i < n; i++) determ *= packmat[i][ld];

  fprintf(fp2,"\n Determinant = ");
  fprintf(fp2,FORMAT_LE, determ);
  
  if (n<=SIZE) {                     
    fprintf(fp2,"\n\n System Solution\n\n");
    for (k=0; k<n; k++) 
      for (i=0; i<n; i++) 
        solution[k]+=matinv[i][k]; 
    WriteVec(fp2, n, solution); 
  }     
                     
  WriteEnd(fp2); 
  fclose(fp2);
               
  printf("\n Results in bandmat5.lst.\n\n");             
  return (0);
}
                                           
// J-P Moreau feb. 1998