#
# Double complex LAPACK routines:  Simple Drivers:  For Linear Equations
# 
# Functions defined here:
#
# ZGESV
# ZGBSV
# ZGTSV
# (lacking ZPOSV, ZPBSV, ZPTSV, ZSYSV, and ZHESV)

#
# Problem:  ZGBSV
#
@PROBLEM zgbsv
@FUNCTION zgbsv
@LANGUAGE FORTRAN
@MAJOR COL
@PATH /LAPACK/Simple/Linear_Equations/
@MOVEABLE 0 
@DESCRIPTION
   ZGBSV computes the solution to a complex system of linear equations
   A * X = B, where A is a band matrix of order N with KL subdiagonals
   and KU superdiagonals, and X and B are N-by-NRHS matrices.
 
   The LU decomposition with partial pivoting and row interchanges is
   used to factor A as A = L * U, where L is a product of permutation
   and unit lower triangular matrices with KL subdiagonals, and U is
   upper triangular with KL+KU superdiagonals.  The factored form of A
   is then used to solve the system of equations A * X = B.
http://www.netlib.org/lapack/
@INPUT 4
@OBJECT SCALAR I kl
KL      (input) INTEGER
   The number of subdiagonals within the band of A.  KL >= 0.
@OBJECT SCALAR I ku
KU      (input) INTEGER
   The number of superdiagonals within the band of A.  KU >= 0.
@OBJECT MATRIX Z ab
AB      (input/output) COMPLEX*16 array, dimension (LDAB,N)
   On entry, the matrix A in band storage, in rows KL+1 to
   2*KL+KU+1; rows 1 to KL of the array need not be set.
   The j-th column of A is stored in the j-th column of the
   array AB as follows:
   AB(KL+KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+KL)
   On exit, details of the factorization: U is stored as an
   upper triangular band matrix with KL+KU superdiagonals in
   rows 1 to KL+KU+1, and the multipliers used during the
   factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
   See below for further details.
@OBJECT MATRIX Z b
B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
   On entry, the N-by-NRHS right hand side matrix B.
   On exit, if INFO = 0, the N-by-NRHS solution matrix X.
@OUTPUT 4
@OBJECT MATRIX Z ab
AB      (input/output) COMPLEX*16 array, dimension (LDAB,N)
 On exit, details of the factorization: U is stored as an
upper triangular band matrix with KL+KU superdiagonals in
rows 1 to KL+KU+1, and the multipliers used during the
factorization are stored in rows KL+KU+2 to 2*KL+KU+1.

@OBJECT VECTOR I ipiv
IPIV    (output) INTEGER array, dimension (N)
   The pivot indices that define the permutation matrix P;
   row i of the matrix was interchanged with row IPIV(i).
@OBJECT MATRIX Z b
B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
   On entry, the N-by-NRHS right hand side matrix B.
   On exit, if INFO = 0, the N-by-NRHS solution matrix X.
@OBJECT SCALAR I info
INFO    (output) INTEGER
   = 0:  successful exit
   < 0:  if INFO = -i, the i-th argument had an illegal value
   > 0:  if INFO = i, U(i,i) is exactly zero.  The factorization
         has been completed, but the factor U is exactly
         singular, and the solution has not been computed.

*  Further Details
*  ===============
*  The band storage scheme is illustrated by the following example, when
*  M = N = 6, KL = 2, KU = 1:
*
*  On entry:                       On exit:
*
*      *    *    *    +    +    +       *    *    *   u14  u25  u36
*      *    *    +    +    +    +       *    *   u13  u24  u35  u46
*      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
*     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
*     a21  a32  a43  a54  a65   *      m21  m32  m43  m54  m65   *
*     a31  a42  a53  a64   *    *      m31  m42  m53  m64   *    *
*
*  Array elements marked * are not used by the routine; elements marked
*  + need not be set on entry, but are required by the routine to store
*  elements of U because of fill-in resulting from the row interchanges.
@COMPLEXITY 2,3
@CALLINGSEQUENCE
#
# n, order of matrix A
@ARG nI2,mI3
# kl, scalar integer
@ARG I0
# ku, scalar integer
@ARG I1
# nrhs, 
@ARG nI3
# AB, matrix of doubles
@ARG I2,O0
# ldab, leading dimension of AB
@ARG lI2
# ipiv, vector of ints, output
@ARG O1
# B, matrix of doubles
@ARG I3,O2
# ldb, leading dimension of B
@ARG lI3
# info, scalar integer
@ARG O3
@COMP mI2=op(+,op(+,array(I1,0),1),op(*,array(I0,0),2))
#
# Code segment
#
@CODE
extern void zgbsv(int *n, int *kl, int *ku, int *nrhs, void *ab,
                  int *ldab, int *ipiv, void *b, int *ldb, int *info);

int ldb,ldab;

#if (*@mI3@ != *@nI2@) return NS_PROT_DIM_MISMATCH;

@O1@ = (int *) malloc(sizeof(int) * (*@nI2@));
*@mO1@ = *@nI2@;

@O3@ = (int *) malloc(sizeof(int));

ldb = (*@mI3@ > 1) ? *@mI3@ : 1;
ldab = (*@mI2@ > 1) ? *@mI2@ : 1;

zgbsv(@nI2@, @I0@, @I1@, @nI3@, @I2@, &ldab, @O1@, @I3@, &ldb, @O3@);

@O0@ = @I2@;
*@nO0@ = *@nI2@;
*@mO0@ = *@mI2@;
@O2@ = @I3@;
*@mO2@ = *@mI3@;
*@nO2@ = *@nI3@;

@END_CODE

#
# ZGTSV
#
@PROBLEM zgtsv
@FUNCTION zgtsv
@LANGUAGE FORTRAN
@MAJOR COL
@PATH /LAPACK/Simple/Linear_Equations/
@MOVEABLE 0 
@DESCRIPTION
   ZGTSV  solves the equation
 
      A*X = B,
 
   where A is an N-by-N tridiagonal matrix, by Gaussian elimination with
   partial pivoting.
 
   Note that the equation  A'*X = B  may be solved by interchanging the
   order of the arguments DU and DL.
http://www.netlib.org/lapack/
@INPUT 4
@OBJECT VECTOR Z dl
DL      (input/output) COMPLEX*16 array, dimension (N-1)
        On entry, DL must contain the (n-1) subdiagonal elements of
        A.
@OBJECT VECTOR Z d
D       (input/output) COMPLEX*16 array, dimension (N)
        On entry, D must contain the diagonal elements of A.
@OBJECT VECTOR Z du
DU      (input/output) COMPLEX*16 array, dimension (N-1)
        On entry, DU must contain the (n-1) superdiagonal elements
        of A.
@OBJECT MATRIX Z b
B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
        On entry, the N-by-NRHS right hand side matrix B.
@OUTPUT 5
@OBJECT VECTOR Z dl
DL      (input/output) COMPLEX*16 array, dimension (N-1)
        On exit, DL is overwritten by the (n-2) elements of the
        second superdiagonal of the upper triangular matrix U from
        the LU factorization of A, in DL(1), ..., DL(n-2).
@OBJECT VECTOR Z d
D       (input/output) COMPLEX*16 array, dimension (N)
        On exit, D is overwritten by the n diagonal elements of U.
@OBJECT VECTOR Z du
DU      (input/output) COMPLEX*16 array, dimension (N-1)
        On exit, DU is overwritten by the (n-1) elements of the first
        superdiagonal of U.
@OBJECT MATRIX Z b
B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
        On exit, if INFO = 0, the N-by-NRHS solution matrix X.
@OBJECT SCALAR I info
INFO    = 0:  successful exit
        < 0:  if INFO = -i, the i-th argument had an illegal value
        > 0:  if INFO = i, U(i,i) is exactly zero, and the solution
              has not been computed.  The factorization has not been
              completed unless i = N.
@COMPLEXITY 2,3
@CALLINGSEQUENCE
# n order of matrix A/vector D.
@ARG mI1, mI3
# nrhs, number cols. of B
@ARG nI3
# dl, array, i/o
@ARG I0,O0
# d, double precision vector, i/o
@ARG I1,O1
# du, double precision vector, i/o
@ARG I2,O2
# b, double precision matrix, i/o
@ARG I3,O3
# ldb, leading dimension of B
@ARG lI3
# info, scalar integer
@ARG O4
@COMP mI0=op(-,mI1,1)
@COMP mI2=op(-,mI1,1)
@CODE
extern void zgtsv(int *n, int *nrhs, void *dl, void *d, void *du,
                  void *b, int *ldb, int *info);

int ldb;

@O4@ = (int *) malloc(sizeof(int));


ldb = (*@mI1@ > 1) ? *@mI1@ : 1;

zgtsv(@mI1@, @nI3@, @I0@, @I1@, @I2@, @I3@, &ldb, @O4@);

@O0@ = @I0@;
*@mO0@ = *@mI0@;

@O1@ = @I1@;
*@mO1@ = *@mI1@;

@O2@ = @I2@;
*@mO2@ = *@mI2@;

@O3@ = @I3@;
*@nO3@ = *@nI3@;
*@mO3@ = *@mI3@;

@END_CODE

#
# ZGESV
#

@PROBLEM zgesv
@FUNCTION zgesv
@LANGUAGE FORTRAN
@MAJOR COL
@PATH /LAPACK/Matrices/LinearSystem/
@MOVEABLE 0 
@DESCRIPTION
From LAPACK -

Compute the solution to a complex system of linear equations
  A * X = b
where A is an N-by-N matrix and X and B are N-by-NRHS matrices.

MATLAB Example : [x y z info ] = netsolve('zgesv',a,b)

http://www.netlib.org/lapack/
@INPUT 2
@OBJECT MATRIX Z a
The N-by-N coefficient matrix A.
@OBJECT MATRIX Z b
The N-by-NRHS matrix of right hand sides B.
@OUTPUT 4
@OBJECT MATRIX Z lu
The factors L and U from the factorization
A = P*L*U; the unit diagonal elements of L are not stored.
@OBJECT VECTOR I ipiv
The pivot indices that define the permutation matrix P;
row i of the matrix was interchanged with row IPIV(i).
@OBJECT MATRIX Z x
The N-by-NRHS solution matrix X.
@OBJECT SCALAR I info
INFO    = 0:  successful exit
        < 0:  if INFO = -i, the i-th argument had an illegal value
        > 0:  if INFO = i, U(i,i) is exactly zero.  The factorization
              has been completed, but the factor U is exactly
              singular, so the solution could not be computed.

@COMPLEXITY 2,3

@CALLINGSEQUENCE
@ARG mI0,nI0,mI1
@ARG nI1
@ARG I0,O0
@ARG lI0
@ARG O1
@ARG I1,O2
@ARG lI1
@ARG O3
@CODE
extern void zgesv(int *n, int *nrhs, void *a, int *lda, int *ipiv,
                  void *b, int *ldb, int *info);

int lda, ldb;

lda =ldb = (*@mI0@ > 1) ? *@mI0@ : 1;

if (*@nI0@ != *@mI0@)
  return NS_PROT_DIM_MISMATCH;
if (*@mI0@ != *@mI1@)
  return NS_PROT_DIM_MISMATCH;

@O1@ = (int *)malloc(sizeof(int)*(*@mI0@));
*@mO1@ = *@mI1@;
@O3@ = (int *)malloc(sizeof(int));

zgesv(@nI0@, @nI1@, @I0@, &lda, @O1@, @I1@, &ldb, @O3@);

@O0@ = @I0@;
*@mO0@ = *@mI0@;
*@nO0@ = *@nI0@;
@O2@ = @I1@;
*@mO2@ = *@mI1@;
*@nO2@ = *@nI1@;

@END_CODE
