/****************************************************************************
 *
 * zcstrd.c
 *
 * Abstract 
 *   Implement the user interface of block Lanczos tridiagonalization
 *   algorithm. It invokes routines implemented in BlkTri.c and LanTri.c 
 *   to perform the tridiagonalization of a complex symmetric matrix.
 *
 * Student:  Guohong Liu
 *
 * ID:       0385117
 *
 ****************************************************************************/

#include <stdio.h>
#include "blklan.h"
 
/***************************************************************************
 *
 * Abstract
 *   Reduce a complex symmetric matrix A stored in normal form  
 *   to complex symmetric tridiagonal form T:
 *      (Q*P)' * A * conjugate(Q*P) = T.   
 *
 * Input  
 *   n      The order of the matrix A.  n >= 0.                
 *   A      Double complex array, dimension (n,n)   
 *   bs     Size of block. b > 0.
 *   S      Double complex array, dimension (n,bs)   
 *          The starting orthonormal columns matrix.
 *           
 * Output  
 *   a      Double complex array, dimension (n)   
 *          The diagonal elements of the tridiagonal matrix T.  
 *   b      Double real array, dimension (n-1)   
 *          The off-diagonal elements of the tridiagonal matrix T.
 *   Q      Double complex array, dimension (n*n)     
 *   P      Double complex array, dimension (n*n)
 *          (Q*P)' * A * conjugate(Q*P) = T   
 *   info   0    Successful exit   
 *          < 0  If info = -i, the i-th argument had an illegal value   
 *          > 0: Exception is thrown
 *
 * Return  
 *   always 0
 *
 ***************************************************************************/

int zcstrd_ (integer *n, doublecomplex *A, integer *bs, doublecomplex *S, 
             doublecomplex *a, doublereal *b, doublecomplex *Q, 
             doublecomplex *P, integer *info)                
{
#ifdef INFO
    DoubleComplexMat M;
#endif

    *info = 0;
    if (!isParamValid (*n, *bs, info)) {
        return 0;
    }
#ifdef INFO
    printf ("\nOrder of matrix A = %d. Block size = %d.\n", *n, *bs);
    M.m = *n;
    M.n = *n;
    M.mat = A;
    showComplexMat (&M, "\nMatrix A");
    M.m = *n;
    M.n = *bs;
    M.mat = S;
    showComplexMat (&M, "\nStarting matrix S");
#endif

    /* Initialize the random number seed.               */
    initRandSeed ();
   
    /* Initialize for block Lanczos tridiagonalization. */
    if (blkTriInit (*n, A, *bs, S, a, b, Q, P)) {
        *info = 1;
        return 0;
    }    

    /* Initialize the Lapack wrapper.                   */
    if (wrapperInit (*n, *bs)) {
        *info = 1;
        return 0;
    }    

    /* Perform block Lanczos tridiagonalization.        */
    blkTri ();
    blkTriEnd ();

    /* Initialize for Lanczos tridiagonalization.       */
    if (lanTriInit ()) {
        *info = 1;
        return 0;
    }
    /* Perform Lanczos tridiagonalization.              */
    lanTri ();
    lanTriEnd ();

    /* Free resources for Lapack wrapper.               */
    freeWrapper ();
    return 0;
}

/***************************************************************************
 *
 * Abstract
 *   Check if parameters given to block Lanczos algorithm are valid.
 *
 * Input
 *   n      The order of the matrix A.  n >= 0.                
 *   bs     Size of block. b > 0.
 *
 * Output     
 *   info   < 0  If info = -i, the i-th argument had an illegal value   
 *          > 0: Exception is thrown
 *
 * Return
 *   TRUE   parameters are valid
 *   FALSE  parameters are invalid
 *
 ***************************************************************************/

Boolean isParamValid (int n, int bs, int *info)
{
    if (n <= 0) {
        *info = -2;
        printf ("\nzcstrd: The order of matrix is less than or equal to 0.\n");
        return FALSE;
    }
    if (bs <= 0) {
        *info = -4;
        printf ("\nzcstrd: The block size is less than or equal to 0.\n");
        return FALSE;
    }    
    if (n % bs != 0) {
        *info = -4;
        printf (
        "\nzcstrd: The order of matrix can not be divided by block size.\n");
        return FALSE;
    }

    return TRUE;
}