/**********************************************************
  *                   MODULE fbando.cpp                     *
  *        Linear banded system without using pivots        *
  * ------------------------------------------------------- *
  * Reference:                                              *
  *                                                         *
  * "Numerical Algorithms with C, By Gisela Engeln-Muellges *
  *  and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].    *
  * ------------------------------------------------------- *
  * See demonstration program tband.cpp                     *
  **********************************************************/
#include "basis.h"

//header of called functions
int banodec(int, int, int, REAL **);
int banosol(int, int, int, REAL **, REAL *);


int bando               /* Linear banded system without using pivots .*/
          (
          int    mod,             /* Modus: 0, 1, 2 ..................*/
          int    n,               /* size of system ..................*/
          int    ld,              /* # of lower co-diagonals .........*/
          int    ud,              /* # of upper co-diagonals .........*/
          REAL ** pmat,           /* condensed input matrix ..........*/
          REAL   *b               /* right hand side .................*/
         )
/*====================================================================*
 *                                                                    *
 *  The function bando solves a linear banded system:  pmat * x = b.  *
 *  Here pmat is a nonsingular n x n matrix in condensed form, i.e.   *
 *  represented in an ld+ud+1 x n matrix for its ld lower and ud upper*
 *  co-diagonals. b denotes the right hand side of the system, and x  *
 *  is the solution.                                                  *
 *                                                                    *
 *  bando uses the Gauss algorithm without column pivot search.       *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Applications:                                                    *
 *   =============                                                    *
 *      Solve linear systems with nonsingular banded system matrices. *
 *      Particularly useful for large sparse and banded and diagonally*
 *      dominant matrices.                                            *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Control parameter:                                               *
 *   ==================                                               *
 *      mod      int mod;                                             *
 *               calling modus for band:                              *
 *       = 0     factor matrix and solve linear system                *
 *       = 1     factor only, store factors in pmat                   *
 *       = 2     solve linear system only; for this call, the factors *
 *               must already be available in pmat, such as when      *
 *               many systems are solved for differing right hand     *
 *               sides and the same system matrix.                    *
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      n        int n;  ( n > 2 )                                    *
 *               Dimension of pmat, size of b                         *
 *      ld       int ld; ( ld >= 0 )                                  *
 *               number of lower co-diagonals                         *
 *      ud       int ud; ( ud >= 0 )                                  *
 *               number of upper co-diagonals                         *
 *      pmat     REAL   *pmat[n];                                     *
 *               mod = 0, 1:                                          *
 *               Matrix of the system in comndensed form.             *
 *               Each row has length at least ld + 1 + ud + min(ld,ud)*
 *               where the columns 0, .., ld-1 denote the lower       *
 *               co-diagonals, column ld stores the diagonal and the  *
 *               columns ld+1, .., ld+ud contain the upper            *
 *               co-diagonals.                                        *
 *               If A is the original uncondensed band matrix, then : *
 *               A[i][k] = pmat[i][ld+k-i],                           *
 *                           for k,i inside the band                  *
 *               mod = 2:                                             *
 *               LU factors in condensed form.                        *
 *      b        REAL   b[n];          ( for mod = 0, 2 )             *
 *               Right hand side                                      *
 *                                                                    *
 *   Output parameters:                                               *
 *   ==================                                               *
 *      pmat     REAL   *pmat[n];      ( for mod = 0, 1 )             *
 *               LU factorization in condensed form                   *
 *      b        REAL   b[n];          ( for mod = 0, 2 )             *
 *               solution vector for the system                       *
 *                                                                    *
 *   Return value :                                                   *
 *   =============                                                    *
 *      = 0      all ok                                               *
 *      = 1      n < 3 or other incorrect input parameter             *
 *      = 2      lack of space                                        *
 *      = 3      Matrix is numerically singular                       *
 *      = 4      wrong calling modus                                  *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Functions used :                                                 *
 *   ================                                                 *
 *      int banodec() : Factor matrix                                 *
 *      int banosol() : solve linear system                           *
 *                                                                    *
 *====================================================================*/
{
  int rc;

  if (n < 3 || ld < 0 || ud < 0 || n < ld + 1 + ud) return (1);

  switch (mod)
  {
    case 0: /* factor matrix and solve system ........................*/
            rc = banodec (n, ld, ud, pmat);
            if (rc == 0)
              return (banosol (n, ld, ud, pmat, b));
            else
              return (rc);

    case 1: /* factor only ...........................................*/
            return (banodec (n, ld, ud, pmat));

    case 2: /* solve only ............................................*/
            return (banosol (n, ld, ud, pmat, b));
  }

  return (3);
}


int banodec             /* Decompose a banded matrix .................*/
            (
             int    n,            /* size of system ..................*/
             int    ld,           /* # of lower co-diagonals .........*/
             int    ud,           /* # of upper co-diagonals .........*/
             REAL * pmat[]        /* condensed input matrix ..........*/
            )
/*====================================================================*
 *                                                                    *
 *  The function banodec factors a condensed banded matrix pmat.      *
 *  banddec uses the Gauss algorithm without column pivot search.     *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      n        int n;  ( n > 2 )                                    *
 *               Dimension of pmat, size of b                         *
 *      ld       int ld; ( ld >= 0 )                                  *
 *               number of lower co-diagonals                         *
 *      ud       int ud; ( ud >= 0 )                                  *
 *               number of upper co-diagonals                         *
 *      pmat     REAL   *pmat[n];                                     *
 *               Matrix of the system in comndensed form.             *
 *               Each row has length at least ld + 1 + ud + min(ld,ud)*
 *               where the columns 0, .., ld-1 denote the lower       *
 *               co-diagonals, column ld stores the diagonal and the  *
 *               columns ld+1, .., ld+ud contain the upper            *
 *               co-diagonals.                                        *
 *                                                                    *
 *   Output parameters:                                               *
 *   ==================                                               *
 *      pmat     REAL   *pmat[n];                                     *
 *               LU factorization in condensed form                   *
 *                                                                    *
 *   Return value :                                                   *
 *   =============                                                    *
 *      = 0      all ok                                               *
 *      = 1      n < 3 or other incorrect input parameter             *
 *      = 2      Matrix is numerically singular                       *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Constants used :     NULL, MACH_EPS                              *
 *   ================                                                 *
 *                                                                    *
 *   Macros: min                                                      *
 *   =======                                                          *
 *                                                                    *
 *====================================================================*/
{
  int      kend, kjend, jm, jk;
  register k, j, i;

  if (ld == 0) return (0);           /* Matrix already in upper       */
                                     /* triangular form               */
  if (pmat == NULL) return (1);

  for (k = 0; k < n; k++)
    if (pmat[k] == NULL) return (1);

  for (i = 0; i < n - 1; i++)               /* loop over all rws      */
  {
    kend =  min (ld + 1, n - i);
    kjend = min (ud + 1, n - i);

    if (ABS(pmat[i][ld]) < MACH_EPS)        /* LU decompsition does   */
      return (2);                           /* not exist              */

    for (k = 1; k != kend; k++)             /* loop over all rows     */
    {                                       /* below row i            */
      pmat[k+i][ld-k] /= pmat[i][ld];

      for (j = 1; j != kjend; j++)
      {
        jk = j + ld - k;
        jm = j + ld;
        pmat[k+i][jk] -= pmat[k+i][ld-k] * pmat[i][jm];
      }
    } /* end k */
  }  /* end i */

  return (0);
}


int banosol             /* Solve a banded system .....................*/
            (
             int    n,            /* size of system ..................*/
             int    ld,           /* # of lower co-diagonals .........*/
             int    ud,           /* # of upper co-diagonals .........*/
             REAL * pmat[],       /* condensed input matrix ..........*/
             REAL   b[]           /* right hand side .................*/
            )
/*====================================================================*
 *                                                                    *
 *  The function banosol solves a factored linear banded system in    *
 *  condensed form using banodec.                                     *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      n        int n;  ( n > 2 )                                    *
 *               Dimension of pmat, size of b                         *
 *      ld       int ld; ( ld >= 0 )                                  *
 *               number of lower co-diagonals                         *
 *      ud       int ud; ( ud >= 0 )                                  *
 *               number of upper co-diagonals                         *
 *      pmat     REAL   *pmat[n];                                     *
 *               Matrices for the factored system in comndensed form. *
 *                                                                    *
 *   Output parameters:                                               *
 *   ==================                                               *
 *      b        REAL   b[n];                                         *
 *               solution vector of linear system                     *
 *                                                                    *
 *   Return value :                                                   *
 *   =============                                                    *
 *      = 0      all ok                                               *
 *      = 1      n < 3 or other incorrect input parameter             *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Constants used :     NULL                                        *
 *   ================                                                 *
 *                                                                    *
 *   Macros: min                                                      *
 *   =======                                                          *
 *                                                                    *
 *====================================================================*/
{
  int      kend;
  register i, k;
                                        /*  Invalid input parameter   */
  if (n < 3 || ld < 0 || ud < 0 || ud + ld + 1 > n) return (1);

  if (pmat == NULL || b == NULL) return (1);

  for (i = 0; i < n; i++)
    if (pmat[i] == NULL) return (1);

  for (i = 0; i < n - 1; i++)
  {
    kend = min (ld + 1, n - i);
    for (k = 1; k != kend; k++)
      b[k+i] -= pmat[k+i][ld-k] * b[i];
  }

  for (i = n - 1; i >= 0 ; i--)          /* back substitution        */
  {
    kend = min (ud + 1, n - i);
    for (k = 1; k < kend; k++)
      b[i] -= pmat[i][k+ld] * b[i+k];
    b[i] /= pmat[i][ld];
  }

  return (0);
}

// -------------------------- END fbando.cpp --------------------------