/* nag_dgbtrs (f07bec) Example Program.

 *

 * Copyright 2001 Numerical Algorithms Group.

 *

 * Mark 7, 2001.

 */



#include <stdio.h>

#include <nag.h>

#include <nag_stdlib.h>

#include <nagf07.h>

#include <nagx04.h>



int main(void)

{

  /* Scalars */

  Integer  i, ipiv_len, j, kl, ku, n, nrhs, pdab, pdb;

  Integer  exit_status=0;

  NagError fail;

  Nag_OrderType order;



  /* Arrays */

  double  *ab=0, *b=0;

  Integer *ipiv=0;



#ifdef NAG_COLUMN_MAJOR

#define AB(I,J) ab[(J-1)*pdab + kl + ku + I - J]

#define B(I,J) b[(J-1)*pdb + I - 1]

  order = Nag_ColMajor;

#else

#define AB(I,J) ab[(I-1)*pdab + kl + J - I]

#define B(I,J) b[(I-1)*pdb + J - 1]

  order = Nag_RowMajor;

#endif



  INIT_FAIL(fail);

  Vprintf("nag_dgbtrs (f07bec) Example Program Results\n\n");



  /* Skip heading in data file */

  Vscanf("%*[^\n] ");

  Vscanf("%ld%ld%ld%ld%*[^\n] ", &n, &nrhs,

         &kl, &ku);

  ipiv_len = n;

#ifdef NAG_COLUMN_MAJOR

  pdab = 2*kl + ku + 1;

  pdb = n;

#else

  pdab = 2*kl + ku + 1;

  pdb = nrhs;

#endif



  /* Allocate memory */

  if ( !(ab = NAG_ALLOC((2*kl+ku+1) * n, double)) ||

       !(b = NAG_ALLOC(nrhs * n, double)) ||

       !(ipiv = NAG_ALLOC(ipiv_len, Integer)) )

    {

      Vprintf("Allocation failure\n");

      exit_status = -1;

      goto END;

    }



  /* Read A from data file */

  for (i = 1; i <= n; ++i)

    {

      for (j = MAX(i-kl,1); j <= MIN(i+ku,n); ++j)

        Vscanf("%lf", &AB(i,j));

    }

  Vscanf("%*[^\n] ");

  /* Read B from data file */

  for (i = 1; i <= n; ++i)

    {

      for (j = 1; j <= nrhs; ++j)

        Vscanf("%lf", &B(i,j));

    }

  Vscanf("%*[^\n] ");



  /* Factorize A */

  /* nag_dgbtrf (f07bdc).

   * LU factorization of real m by n band matrix

   */

  nag_dgbtrf(order, n, n, kl, ku, ab, pdab, ipiv, &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from nag_dgbtrf (f07bdc).\n%s\n", fail.message);

      exit_status = 1;

      goto END;

    }

  /* Compute solution */

  /* nag_dgbtrs (f07bec).

   * Solution of real band system of linear equations,

   * multiple right-hand sides, matrix already factorized by

   * nag_dgbtrf (f07bdc)

   */

  nag_dgbtrs(order, Nag_NoTrans, n, kl, ku, nrhs, ab, pdab, ipiv,

             b, pdb, &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from nag_dgbtrs (f07bec).\n%s\n", fail.message);

      exit_status = 1;

      goto END;

    }

  /* Print solution */

  /* nag_gen_real_mat_print (x04cac).

   * Print real general matrix (easy-to-use)

   */

  nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n, nrhs, b,

                         pdb, "Solution(s)", 0, &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from nag_gen_real_mat_print (x04cac).\n%s\n",

              fail.message);

      exit_status = 1;

      goto END;

    }

 END:

  if (ab) NAG_FREE(ab);

  if (b) NAG_FREE(b);

  if (ipiv) NAG_FREE(ipiv);

  return exit_status;

}